เตรียมตัวก่อนเรียน

บทเรียนนี้จะกล่าวถึงการวิเคราะห์โมเดลพหุระดับด้วยวิธีการแบบเบส์ ซึ่งมีความแตกต่างไปจากการวิเคราะห์ด้วยวิธีการแบบดั้งเดิมที่เรียกว่าวิธีการแบบ frequentist ส่วนแรกของบทเรียนจะปูพื้นฐานเกี่ยวกับสถิติแบบเบส์ก่อน จากนั้นจึงกล่าวถึงการใช้สถิติแบบเบส์ในการวิเคราะห์โมเดลพหุระดับ

โปรแกรมที่ใช้ในการบรรยายประกอบด้วย

  • Jamovi

  • R

  • JAGS (Just Another Gibb Sampler)

  • STAN

ทั้งนี้เพื่อให้ไม่เป็นการเสียเวลา ขอให้ผู้เรียนดาวน์โหลดและติดตั้งโปรแกรมดังกล่าวให้พร้อมก่อนวันที่เรียน

ความรู้พื้นฐาน

หัวข้อนี้เป็นความรู้พื้นฐานเกี่ยวกับสถิติแบบดั้งเดิม และ สถิติแบบเบส์ โดยจะกล่าวถึงข้อตกลงเบื้องต้นของสถิติแต่ละประเภท และข้อจำกัดของสถิติแบบดั้งเดิมที่สามารถใช้สถิติแบบเบส์เข้ามาช่วยแก้ปัญหาได้ รายละเอียดมีดังนี้

สถิติแบบ Frequentist

สถิติแบบดั้งเดิมหรือที่เรียกกันว่า frequentist เป็นสถิติวิเคราะห์ที่ถูกใช้งานเป็นหลักในปัจจุบัน สถิติเกือบทุกตัวที่ผู้เรียนได้เรียนมาจนถึงปัจจุบันนี้ล้วนเป็นสถิติวิเคราะห์ภายใต้กระบวนทัศน์แบบ frequentist ทั้งสิ้น สถิติแบบ frequentist นี้มีข้อตกลงเบื้องต้นว่า (1) พารามิเตอร์เป็นค่าคงที่ ที่มีอยู่จริงในประชากรเป้าหมาย อย่างไรก็ตามด้วยข้อจำกัดที่ไม่สามารถเก็บรวบรวมข้อมูลจากทั้งประชากรได้ทำให้ผู้วิเคราะห์ไม่ทราบค่าพารามิเตอร์ดังกล่าว และจะประมาณโดยใช้ข้อมูลจากตัวอย่าง (2) ข้อมูลตัวอย่าง (sample data) เขียนแทนด้วย \(\bf{x}= \left \{\bf{x_1}, \bf{x_2},...,\bf{x_n} \right \}\) เป็นค่าที่สุ่มมาจากการประชากร ดังนั้นตัวอย่างที่สุ่มได้จึงเป็นค่าสุ่มที่และมีความเป็นไปได้ทั้งหมดเท่ากับ เมื่อ \(N\) คือขนาดประชากร และ \(n\) คือขนาดตัวอย่าง จากข้อตกลงเบื้องต้นข้อที่ (2) นี้จึงทำให้ (3) ค่าสถิติที่ประมาณได้จากข้อมูลตัวอย่างเป็นตัวแปรสุ่มที่มีการแจกแจงความน่าจะเป็น เรียกการแจกแจงความน่าจะเป็นของค่าสถิตินี้ว่า การแจกแจงความน่าจะเป็นของฟังก์ชันที่ได้จากตัวอย่างสุ่ม (sampling distribution)

ลักษณะการแจกแจงความน่าจะเป็นของค่าสถิติดังกล่าวมักประมาณได้ด้วย 2 วิธีการ วิธีการแรก ประมาณโดยใช้ทฤษฎีลิมิตลู่เข้าสู่ส่วนกลาง (central limit theorem) ซึ่งเป็นการคาดการณ์ลักษณะการแจกแจงของค่าสถิติในเชิงทฤษฎี เช่น การประมาณช่วงความเชื่อมั่นของค่าเฉลี่ยมีสูตรเป็น \(\overline{X}\pm t_{\alpha/2,n-1} \frac{S}{\sqrt{n}}\) จะเห็นว่าการประมาณช่วงความเชื่อมั่นดังกล่าวมีการใช้การแจกแจงที (student t distribution) ในการคำนวณส่วน margin of error การแจกแจงดังกล่าวพิสูจน์ในทางทฤษฎีโดยใช้ central limit theorem ได้ว่า เป็น sampling distribution ของ \(\overline{X}\) หากข้อมูลตัวอย่างที่ใช้ในการวิเคราะห์สุ่มมาจากประชากรที่มีการแจกแจงแบบปกติหรือใกล้เคียง แต่ไม่ทราบค่าความแปรปรวนของประชากร หรือถ้าหากประชากรไม่ได้มีการแจกแจงแบบปกติ แต่ตัวอย่างมีขนาดใหญ่เพียงพอ ก็ยังสามารถใช้ sampling distribution ดังกล่าวประมาณการแจกแจงของ \(\overline{X}\) ได้อยู่ หรือในการทดสอบสมมุติฐาน \(H_0: \mu=\mu_0\) สถิติทดสอบที่ใช้คือ \(t^*=\frac{\overline{X}-\mu_0}{S/ \sqrt{n}}\) จะมี sampling distributionในทางทฤษฎีเป็นการแจกแจงแบบที ที่มีองศาความเป็นอิสระเท่ากับ \(n-1\) เช่นเดียวกันหากข้อกำหนดเบื้องต้นของ central limit theorem เป็นจริงหรือใกล้เคียง การทราบการแจกแจงของค่าสถิติดังกล่าวทำให้ผู้วิเคราะห์สามารถประเมินความคลาดเคลื่อนในการตัดสินใจของการทดสอบสมมุติฐาน และสรุปผลได้ดังที่ได้เรียนไปแล้วในรายวิชาพื้นฐาน

ความถูกต้องของทฤษฎีนี้ขึ้นอยู่กับปัจจัยแวดล้อมที่เกี่ยวข้องว่ามีความสอดคล้องกับทฤษฎีมากน้อยเพียงใด ปัจจัยหนึ่งที่มีความสำคัญคือขนาดตัวอย่าง กล่าวคือหากขนาดตัวอย่างมีค่าน้อยเกินไป หรือการแจกแจงของข้อมูลมีความเบี่ยงเบนออกไปจากข้อตกลงเบื้องต้นของการวิเคราะห์ไปมาก ค่าสถิติที่คำนวณได้อาจไม่มีลักษณะการแจกแจงที่เป็นไปตามทฤษฎี และส่งผลให้การอนุมานเชิงสถิติมีความคลาดเคลื่อน อีกวิธีการหนึ่งที่สามารถใช้ในการประมาณ sampling distribution ของค่าสถิติได้คือใช้วิธีการสุ่มซ้ำ (resampling method) เช่น boostraping หรือ jackknife เพื่อประมาณแนวโน้มการแจกแจงของค่าสถิติ เป็นต้น


t.test(dat$Ach, mu=5.5)
## 
##  One Sample t-test
## 
## data:  dat$Ach
## t = 5.078, df = 149, p-value = 1.123e-06
## alternative hypothesis: true mean is not equal to 5.5
## 95 percent confidence interval:
##  5.709732 5.976934
## sample estimates:
## mean of x 
##  5.843333

คำถาม

  1. output ข้างต้นแสดงผลการประมาณช่วงความเชื่อมั่น และทดสอบสมมุติฐานทางสถิติเกี่ยวกับผลสัมฤทธิ์ทางการเรียนของนักเรียน (คะแนนเต็ม 10) โดยจากการประมาณด้วยช่วงความเชื่อมั่น 95% ในข้างต้นพบว่ามีค่าเท่ากับ \([5.71, 5.98]\) ซึ่งหมายความว่าอะไร?

  2. จากผลการทดสอบสมมุติฐาน \(H_0: \mu_{Ach}=5.5\) พบว่ามีค่า p-value < .000 หมายความว่าอะไร?

  3. เราสามารถตัดสินใจยอมรับสมมุติฐาน \(H_0\) ได้หรือไม่?


จากข้อตกลงเบื้องต้นของสถิติแบบ frequentist นี้ทำให้การวิเคราะห์ข้อมูลเกิดข้อจำกัดต่าง ๆ ดังนี้

  1. ด้านที่สำคัญคือด้านการแปลผลการวิเคราะห์ที่จะเห็นว่าความหมายของผลการวิเคราะห์ต่าง ๆ เป็นการกล่าวถึงพารามิเตอร์ที่สนใจแบบอ้อม ๆ ทั้งสิ้น ไม่มีการวิเคราะห์ใดที่สามารถอ้างอิงไปยังพารามิเตอร์ที่สนใจได้โดยตรง ทั้งนี้เป็นเพราะการอนุมานเชิงสถิติแบบ frequentist พึ่งพาเครื่องมือหลักคือ samping distribution ของค่าสถิติ ในขณะที่พารามิเตอร์ถูกสมมุติให้เป็นค่าคงที่ตั้งแต่ต้น probability statement ต่าง ๆ จึงเป็นของค่าสถิติทั้งสิ้น ดังความความเชื่อมั่นที่รับประกันในช่วงความเชือมั่น และค่าความคลาดเคลื่อนในการทดสอบสมมุติฐานจึงเป็นค่าความน่าจะเป็นที่ใช้รับประกันความเป็นไปได้ของค่าสถิติทั้งสิ้นไม่ใช่ของค่าพารามิเตอร์ที่สนใจจริง ๆ ในการวิเคราะห์

  2. ด้านที่สองคือด้านการประมาณค่าพารามิเตอร์ภายใต้โมเดลซับซ้อน เครื่องมือสำคัญสำหรับประมาณค่าพารามิเตอร์ในโมเดลต่าง ๆ ของสถิติแบบ frequentist คือตัวประมาณ เช่น ตัวประมาณกำลังสองน้อยสุด ตัวประมาณภาวะความควรจะเป็นสูงสุด ซึ่งล้วนเป็นวิธีการที่จะหาค่าประมาณที่ดีที่สุดโดยอิงจากจุดต่ำสุด หรือจุดสูงสุดของฟังก์ชันวัตถุประสงค์ ดังตัวอย่างในรูปด้านล่าง ซึ่งในกรณีที่โมเดลมีความซับซ้อนมาก ๆ ฟังก์ชันวัตถุประสงค์ดังกล่าวจะมีความซับซ้อนโดยมีจำนวนมิติที่สูงขึ้น และอาจมีจุดสูงต่ำ และ/หรือ จุดต่ำสุดสัมพัทธ์จำนวนมาก ซึ่งเป็นอุปสรรคให้การประมาณค่าด้วยวิธีการในข้างต้นอาจได้ค่าประมาณที่คลาดเคลื่อนไปจากค่าที่เหมาะสมที่สุด หรืออัลกอริทึมไม่สามารถหาค่าประมาณที่เหมาะสมได้ภายใต้จำนวนรอบของการประมาณที่กำหนด


  1. สืบเนื่องจากที่การอนุมานเชิงสถิติแบบ frequentist นั้นมักอิงกับทฤษฎีความน่าจะเป็นที่อิงข้อกำหนดเบื้องต้นที่ขนาดตัวอย่างต้องมีขนาดใหญ่เพียงพอ ซึ่งทำให้สถิติแบบ frequentist มักทำงานได้ไม่ดีในกรณีที่ตัวอย่างมีขนาดเล็ก

สถิติแบบเบส์ (Bayesian Statistics)

สถิติแบบเบส์เป็นกระบวนทัศน์ในการวิเคราะห์ข้อมูลอีกลักษณะหนึ่งที่มีความแตกต่างจากไปกระบวนทัศน์แบบดั้งเดิมที่นักวิเคราะห์ข้อมูลทั่วไปคุ้นเคยกัน ซึ่งทำให้สถิติแบบเบส์มีจุดเด่นและสามารถแก้ไขข้อจำกัดที่พบในสถิติแบบ frequentist ได้ ความแตกต่างดังกล่าวเริ่มตั้งแต่ข้อตกลงเบื้องต้นของสถิติแบบเบส์ที่แตกต่างจากสถิติแบบดั้งเดิมอย่างมาก กล่าวคือ

(1) พารามิเตอร์เป็นตัวแปรสุ่ม สถิติแบบเบส์มองว่าพารามิเตอร์เป็นค่าที่ผู้วิเคราะห์สนใจจะศึกษาและต้องการทราบอย่างแท้จริง แต่อย่างไรก็ตามเนื่องจากผู้วิเคราะห์ไม่ทราบข้อมูลประชากร การคาดการณ์หรือความเชื่อ (belief) เกี่ยวกับพารามิเตอร์ดังกล่าวจึงมีความไม่แน่นอน สภาพดังกล่าวจึงสมเหตุสมผลที่จะกำหนด probability statement ให้กับพารามิเตอร์ในโมเดล

จากข้อกำหนดนี้ผลการประมาณค่าพารามิเตอร์ที่ได้จากวิธีการแบบเบส์จึงจะไม่ได้ค่าประมาณแบบค่าเดียว แต่จะได้เป็นการประมาณการแจกแจงของพารามิเตอร์แทน เช่น การประมาณค่าเฉลี่ยผลสัมฤทธิ์ของนักเรียน ผลการวิเคราะห์แบบเบส์จะไม่ได้ให้เป็นค่าประมาณเพียงค่าเดียวคือ 5.84 คะแนน ดัง output ข้างต้น แต่จะให้เป็นการแจกแจงของพารามิเตอร์ค่าเฉลี่ยผลสัมฤทธิ์ทางการเรียน ซึ่งเรียกว่าการแจกแจงความน่าจะเป็นภายหลัง (posterior distribution) และสารสนเทศเกี่ยวกับค่าเฉลี่ยของผลสัมฤทธิ์ดังกล่าวจะถูกสกัดออกมาจากการแจกแจงนี้

(2) ข้อมูลตัวอย่างเป็นค่าคงที่ ถึงแม้ว่าในความเป็นจริงตัวอย่างที่สุ่มมาจากประชากรจะมีความเป็นไปได้ที่หลากหลาย และมีความไม่แน่นอนในการได้มาซึ่งตัวอย่างแต่ละชุด อย่างไรก็ตามเมื่อดำเนินการสุ่มตัวอย่างมาทำการวิเคราะห์แล้ว ความเป็นไปได้จำนวนมากก่อนการสุ่มตัวอย่างจะเหลือเพียงความเป็นไปได้เดียว สถิติแบบเบส์จึงถือว่าข้อมูลตัวอย่างเป็นค่าคงที่ ไม่มีการกำหนด probability statetment ให้กับข้อมูลตัวอย่างนี้

สถิติแบบเบส์สามารถนำมาใช้เพื่อแก้ไขข้อจำกัดที่เกิดขึ้นจากสถิติแบบ frequentist ได้ดังนี้

  1. จากข้อกำหนดเบื้องต้นข้างต้น การอนุมานเชิงสถิติด้วยสถิติแบบเบส์จะดำเนินการผ่านการแจกแจงความน่าจะเป็นของพารามิเตอร์โดยตรงที่เรียกว่า การแจกแจงความน่าจะเป็นภายหลัง (posterior distribution) การแจกแจงดังกล่าวเป็นเครื่องมือสำคัญในการอนุมานเชิงสถิติแบบเบส์ แทน sampling distribution ในการอนุมานเชิงสถิติแบบดั้งเดิม ซึ่งจุดนี้ทำให้สถิติแบบเบส์สร้างความแตกต่างอย่างมากเมื่อเปรียบเทียบกันสถิติแบบดั้งเดิม โดยจะเห็นว่าการประมาณค่าพารามิเตอร์แบบเบส์ไม่ได้อิงกับค่าสถิติเพียงค่าเดียวเหมือนกับสถิติแบบดั้งเดิม แต่อิงกับการแจกแจงความน่าจะเป็นภายหลัง

  2. การที่สถิติแบบเบส์ใช้การแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์เป็นเครื่องมือหลักในการอนุมานเชิงสถิติ จึงทำให้การวิเคราะห์ทำได้อย่างตรงไปตรงมาโดยไม่ต้องอาศัยทฤษฎีหรือวิธีการที่ซับซ้อนเหมือนกับสถิติแบบดั้งเดิม กล่าวคือผู้วิเคราะห์ใช้เพียงสถิติพื้นฐาน เช่น ค่าเฉลี่ย มัธยฐาน ฐานนิยม ส่วนเบี่ยงเบนมาตรฐาน หรือเปอร์เซ็นไทล์ เป็นเครื่องมือในการสรุปสารสนเทศเกี่ยวกับพารามิเตอร์ที่สนใจออกมาจากการแจกแจงความน่าจะเป็นภายหลังที่ประมาณได้ ก็เพียงพอแล้ว นอกจากนี้การแปลผลเพื่ออนุมานเกี่ยวกับพารามิเตอร์ยังทำได้ง่ายและตรงไปตรงมา และได้ข้อสรุปโดยตรงไปยังปริภูมิของพารามิเตอร์ที่เป็นเป้าหมาย แตกต่างจากข้อสรุปที่ได้จากสถิติแบบดั้งเดิมที่เป็นข้อสรุปเกี่ยวกับการแจกแจงของค่าสถิติซึ่งอยู่บนปริภูมิตัวอย่างเท่านั้น

  3. การประมาณการแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ในโมเดล จะใช้สารสนเทศสองส่วนด้วยกัน ส่วนแรกเป็นสมมุติฐานหรือความเชื่อเบื้องต้น (prior belief) เกี่ยวกับพารามิเตอร์ที่สนใจ ซึ่งกำหนดอยู่ในรูปการแจกแจงความน่าจะเป็นที่เรียกว่า การแจกแจงความน่าจะเป็นก่อนหน้าของพารามิเตอร์ (prior distribution) และส่วนที่สองคือสารสนเทศที่ได้จากข้อมูลตัวอย่าง ซึ่งแตกต่างจากสถิติแบบดั้งเดิมที่จะใช้สารสนเทศจากข้อมูลตัวอย่างเท่านั้น

ลองพิจารณาตัวอย่างต่อไปนี้

สมมุติว่าต้องการคาดการณ์คะแนนสอบของนายบูม ที่มีความเป็นไปได้ 4 แบบดังรูป

ก่อนการเก็บรวบรวมข้อมูลจริงในชั้นเรียนของนายบูม สมมุติว่าผู้วิเคราะห์มีข้อมูลในอดีตที่เกี่ยวข้องกับคะแนนสอบของนายบูมดังนี้

  1. ไม่มีข้อมูลใด ๆ ที่เป็นประโยชน์ (ให้น้ำหนักกับความเป็นไปได้ A1, A2, A3 และ A4 อย่างไร?)

  2. หากทราบเพิ่มเติมว่าการสอบที่ผ่าน ๆ มาโดยเฉลี่ยนายบูมสอบได้คะแนนคิดเป็นอันดับที่ 4 ของชั้นเรียนนี้ (น้ำหนักของความเป็นไปได้ทั้ง 4 จะเหมือนเดิมมั้ย?)

  3. หากทราบเพิ่มเติมอีกว่า การสอบที่ผ่าน ๆ มา แก้วซึ่งเป็นเพื่อนในห้องเรียนของนายบูม สอบได้คะแนน 38 คะแนน คิดเป็นคะแนนที่น้อยที่สุดในชั้นเรียนนี้ (น้ำหนักของความเป็นไปได้ทั้ง 4 เหมือนเดิมหรือไม่?)

  4. หากทราบอีกว่า ในการสอบที่ผ่าน ๆ มา นิดสอบได้คะแนนโดยเฉลี่ยคิดเป็น 89 คะแนน ซึ่งเป็นอันดับที่ 3 ของชั้นเรียนนี้ (น้ำหนักของความเป็นไปได้ทั้ง 4 เหมือนเดิมหรือไม่?)

การแจกแจงก่อนหน้าของคะแนนสอบนายบูม

ผู้วิเคราะห์สามารถปรับสารสนเทศในการแจกแจงความน่าจะเป็นก่อนหน้าในข้างต้นได้ หากมีการเก็บรวบรวมข้อมูลตัวอย่างเพิ่มเติม เช่น

สมมุติว่าการแจกแจงความน่าจะเป็นก่อนหน้าของคะแนนนายบูมเป็นแบบ uniform และมีข้อมูลการสอบครั้งล่าสุดเพิ่มเติมดังรูป

อีกตัวอย่างหนึ่ง สมมุติว่า การแจกแจงความน่าจะเป็นก่อนหน้าของคะแนนนายบูมมีน้ำหนักสูงที่ A3 และมีข้อมูลการสอบครั้งล่าสุดเพิ่มเติมดังรูป

ทฤษฎีบทของเบส์ (Bayes’ theorem)

คำถามที่เกิดขึ้นจากตัวอย่างในข้างต้นคือ จะสามารถจัดสรรหรือปรับปรุงน้ำหนักให้กับแต่ละความเป็นไปได้ในการแจกแจงความน่าจะเป็นก่อนหน้าอย่างไร เมื่อมีข้อมูลเพิ่มเติม คำตอบคือ เราสามารถใช้ทฤษฎีบทของเบส์เพื่อดำเนินการดังกล่าวได้ ทฤษฎีบทของเบส์มีรากฐานมาจากการคำนวณความน่าจะเป็นแบบมีเงื่อนไขดังนี้ \(P(A|B)=P(A \cap B)/P(B)\)

ดังนั้นหากกำหนดให้ \(\theta\) คือพารามิเตอร์ภายในโมเดล และ \(\bf{x}\) คือข้อมูลตัวอย่าง การแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ที่ต้องการคือ \(p(\theta|\bf{x})\) ซึ่งสามารถคำนวณได้จาก

\(p(\theta|\bf{x})=\frac{p(\theta,\bf{x})}{p(\bf{x})}\)

จากสมการในรูปข้างต้น เราสามารถปรับรูปสมการใหม่ได้ดังนี้

\(p(\theta|\bf{x})p(\bf{x})=p(\theta,\bf{x})\) ——- (1)

ในทางกลับกันจากสมการความน่าจะเป็นแบบมีเงื่อนไขจะได้ว่า \(p(\bf{x}|\theta)=p(\theta, \bf{x})/p(\bf{\theta})\) ดังนั้น

\(p(\bf{x}|\theta)p(\bf{\theta})=p(\theta,\bf{x})\) ——- (2)

จากสมการที่ (1) กับ (2) จะได้ว่า

\(p(\theta|\bf{x})p(\bf{x})=p(\bf{x}|\theta)p(\bf{\theta})\)

ดังนั้นจะได้ว่าการแจกแจงความน่าจะเป็นภายหลังสามารถคำนวณได้จาก

\(p(\bf{\theta}|\bf{x})=\frac{p(\bf{x}|\theta)p(\theta)}{p(\bf{x})}\)

เมื่อ \(p(\bf{\theta}|\bf{x})\) คือการแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ \(\theta\) จะเห็นว่าเป็นการแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ดังกล่าวเมื่อกำหนดข้อมูลตัวอย่าง \(\bf{x}\) ส่วน \(p(\bf{x}|\theta)\) คือฟังก์ชันภาวะความควรจะเป็น \(p(\theta)\) คือฟังก์ชันความน่าจะเป็นก่อนหน้าของพารามิเตอร์ (prior distribution) และ \(p(\bf{x})\) คือค่าคงที่

สมการข้างต้นจึงสามารถเขียนใหม่ได้ดังนี้

\(p(\bf{\theta}|\bf{x}) \propto p(\bf{x}|\theta)p(\theta)\)

จากสมการในข้างต้นจะเห็นว่าการแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ เกิดจากการนำสารสนเทศสองส่วนเข้ามาประมวลผลร่วมกัน ได้แก่ สมมุติฐานหรือความเชื่อเบื้องต้น (prior belief) เกี่ยวกับพารามิเตอร์ในโมเดล ที่กำหนดผ่านการแจกแจงความน่าจะเป็นก่อนหน้า และสารสนเทศจากข้อมูลตัวอย่างภายในฟังก์ชันภาวะความควรจะเป็น

ตัวอย่างการอนุมานเชิงสถิติแบบเบส์

โดยปกติขั้นตอนสำคัญของการอนุมานเชิงสถิติแบบเบส์มี 5 ขั้นตอน ได้แก่

  1. การระบุขอบเขตด้านตัวแปร

  2. กำหนดโมเดลของค่าสังเกต

  3. กำหนดการแจกแจงความน่าจะเป็นเบื้องต้นของพารามิเตอร์

  4. ประมาณการแจกแจงความน่าจะเป็นภายหลัง

  5. วิเคราะห์ผลลัพธ์ที่ได้จากการแจกแจงความน่าจะเป็นภายหลัง

หัวข้อนี้จะแสดงตัวอย่างการอนุมานเชิงสถิติแบบเบส์ โดยใช้ตัวอย่างง่าย ๆ ดังนี้

กำหนดตัวแปรในการวิเคราะห์

สมมุติว่าผู้วิจัยต้องการวิเคราะห์ความลำเอียงในการออกหน้าหัวหรือก้อยของเหรียญที่ผลิตขึ้นมารุ่นหนึ่ง การวิเคราะห์นี้ตัวแปรที่สนใจจึงเป็นผลลัพธ์ที่สังเกตได้จากการโยนเหรียญในแต่ละครั้ง กำหนดให้ \(y_i\) แทนค่าสังเกตที่ได้จากการโยนเหรียญในครั้งที่ \(i\) จะได้ว่า

\(y_i=\left\{\begin{matrix} 1 & ; H\\ 0 & ; T \end{matrix}\right.\)

กำหนดโมเดลของค่าสังเกต

เมื่อกำหนดตัวแปรที่สนใจแล้วจะพบว่าตัวแปรที่สนใจมีเพียงตัวแปรตาม 1 ตัว โดยมีค่าสังเกตที่เป็นไปได้เพียง 2 ค่าได้แก่ H หรือ T ซึ่งแทนด้วย 1 กับ 0 ตามลำดับ โมเดลของค่าสังเกต \(y_i\) นี้จึงสามารถกำหนดได้ด้วยการแจกแจงความน่าจะเป็นแบบ Bernoulli กล่าวคือ ความน่าจะเป็นที่จะออกหน้าหัวมีค่าเท่ากับ \(p(y_i=1|\theta)=\theta\) และความน่าจะเป็นที่ออกหน้าก้อยมีค่าเท่ากับ \(p(y_i=0|\theta)=1-\theta\) ซึ่งเขียนเป็นสมการรวมได้ดังนี้

\(p(y_i|\theta)=\theta^{y_i}(1-\theta)^{1-y_i}\) เมื่อ \(i = 1,2,3,...,n\)

จากโมเดลในข้างต้นจะเห็นว่าการเกิดค่าสังเกต \(y_i\) ขึ้นอยู่กับพารามิเตอร์ \(\theta\) ที่มีค่าอยู่บนช่วง \([0,1]\) หากพารามิเตอร์ดังกล่าวมีค่าเท่ากับ \(\theta=0.5\) นั่นหมายถึงโอกาสในการออกหน้าหัวและก้อยมีค่าเท่ากัน แต่ถ้าหากพารามิเตอร์ \(\theta>0.5\) นั่นหมายถึงโอกาสในการออกหน้าหัวมีมากกว่า จากความหมายดังกล่าวทำให้สามารถมองได้ว่าพารามิเตอร์ \(\theta\) เป็นพารามิเตอร์ความลำเอียงของเหรียญ หากประมาณพารามิเตอร์ดังกล่าวได้จะสามารถสรุปความลำเอียงของเหรียญได้

เนื่องจากเราไม่ได้เก็บรวบรวมข้อมูลเพียงค่าสังเกตเดียว แต่เราทำการทดลองซ้ำ ๆ จำนวน \(n\) ครั้ง เนื่องจากการทดลองแต่ละครั้งเป็นอิสระซึ่งกันและกัน ทำให้โมเดลของชุดค่าสังเกตหรือเรียกอย่างเป็นทางการว่า ฟังก์ชันภาวะความควรจะเป็นดังนี้

\(p(\bf{y}|\theta)=p(y_1|\theta)p(y_2|\theta)...p(y_n|\theta)\)

ซึ่งมีค่าเท่ากับ

\(p(\bf{y}|\theta)=\theta^{\sum_{i=1}^ny_i}(1-\theta)^{n-\sum_{i=1}^ny_i}\)

ในฟังก์ชันภาวะความควรจะเป็นนี้จะให้ค่าความน่าจะเป็นหรือค่าความหนาแน่นของข้อมูลตัวอย่าง \(\bf{y}\) บนแต่ละค่าที่เป็นไปได้ของพารามิเตอร์ \(\theta\) เนื่องจากข้อมูลตัวอย่างเป็นค่าคงที่ ดังนั้นค่าความหนาแน่นนี้จึงเปลี่ยนแปลงไปตามค่าของพารามิเตอร์ \(\theta\) โดยหากค่าความหนาแน่นนี้มีค่าสูง ณ ค่า \(\theta\) ค่าใด นั่นหมายถึงว่าพารามิเตอร์ค่านั้นทำให้โมเดลของค่าสังเกตกับค่าสังเกตจริงมีความสอดคล้องกันสูง ข้อสังเกตนี้ถูกนำไปใช้เพื่อหาค่าประมาณด้วยวิธีภาวะความควรจะเป็นสูงสุด (maximum likelihood: ML) แต่ในวิธีการแบบเบส์จะนำสารสนเทศส่วนนี้ไปประมวลร่วมกับความเป็นไปได้ของพารามิเตอร์ในการแจกแจงความน่าจะเป็นเบื้องต้นก่อน

กำหนดการแจกแจงความน่าจะเป็นเบื้องต้นของพารามิเตอร์

การกำหนดการแจกแจงความน่าจะเป็นเบื้องต้นของพารามิเตอร์ \(\theta\) สามารถทำได้หลายลักษณะ ทั้งแบบไม่ต่อเนื่อง และแบบต่อเนื่อง ขึ้นอยู่กับขอบเขตของการวิเคราะห์ ข้อมูลในอดีตที่ใช้สนับสนุน และธรรมชาติของพารามิเตอร์นั้น ในตัวอย่างข้างต้นพารามิเตอร์ \(\theta\) ต้องมีค่าอยู่บนช่วง \([0,1]\)

เพื่อให้ตัวอย่างนี้ทำความเข้าใจได้ง่าย จึงเลือกกำหนดการแจกแจงความน่าจะเป็นเบื้องต้นของพารามิเตอร์ \(\theta\) เป็นการแจกแจงแบบไม่ต่อเนื่อง โดยมีค่าที่เป็นไปได้เท่ากับ 0.0, 0.1, 0.2, 0.3, …., 1.0 และมีการแจกแจงของความน่าจะเป็นก่อนหน้าดังรูป

# prior distribution
theta<-seq(0,1,0.1)
p.theta<-c(0.01,0.04,0.075,0.1,0.15,0.25,0.15,0.1,0.075,0.04,0.01)
plot(theta,p.theta, 
     xlab=expression(theta), 
    ylab=expression(p(theta)), 
    type="h", lwd=2, col="orange",
    main="Prior Distribution")
points(theta,p.theta, type="p",pch=16, cex=2, col="orange")

การประมาณการแจกแจงความน่าจะเป็นภายหลัง

จากการทดลองโยนเหรียญจำนวน 10 ครั้งพบว่า ออกหน้าหัวจำนวน 6 ครั้ง ดังนั้นจะได้ว่า

\(p(\bf{y}|\theta)=\theta^{6}(1-\theta)^4\)

จากขอบเขตของ \(theta\) ที่กำหนดในการแจกแจงความน่าจะเป็นเบื้องต้นจะได้ว่าฟังก์ชันภาวะความควรจะเป็น มีลักษณะดังรูป

# likelihood function
L<-theta^6*(1-theta)^4
plot(theta,L, 
     xlab=expression(theta), 
    ylab="p(y|theta)",
    type="h", lwd=2, col="orange",
    main="Likelihood Function")
points(theta,L, type="p",pch=16, cex=2, col="orange")

จากทฤษฎีบทของเบส์เราสามารถรวมสารสนเทศจากการแจกแจงความน่าจะเป็นเบื้องต้น กับฟังก์ชันภาวะความควรจะเป็นได้ดังนี้

# marginal y
p.y<-sum(p.theta*L)

par(mfrow=c(1,2))
# posterior of theta given y
post.theta<-L*p.theta/p.y
plot(theta,post.theta, 
     xlab=expression(theta), 
    ylab="p(theta|y)",
    type="h", lwd=2, col="blue",
    ylim=c(0,0.5),
    main="Posterior Distribution")
points(theta,post.theta, type="p",pch=16, cex=2, col="blue")

# prior and posterior in the same plot
plot(theta,p.theta, 
     xlab=expression(theta), 
    ylab="Probability", 
    type="l", lwd=2, col="orange",
    ylim=c(0,0.5))
points(theta,p.theta, type="p",pch=16, cex=2, col="orange")

points(theta,post.theta, 
    type="l", lwd=2, col="blue")
points(theta,post.theta, type="p",pch=16, cex=2, col="blue")
legend(-0.05,0.48, legend=c("Prior","Posterior"), col=c("orange","blue"), lty=1, bty="n", lwd=2)

การอนุมานเชิงสถิติจากการแจกแจงความน่าจะเป็นภายหลัง

เมื่อได้การแจกแจงความน่าจะเป็นภายหลังที่ต้องการแล้ว ผู้วิเคราะห์สามารถวิเคราะห์การแจกแจงดังกล่าวเพื่ออนุมานเกี่ยวกับพารามิเตอร์ความลำเอียงที่ต้องการศึกษาได้ โดยสามารถทำได้หลายลักษณะ ทั้งการประมาณค่าแบบจุด การประมาณค่าแบบช่วง และการทดสอบสมมุติฐาน รายละเอียดดังนี้

(1) การประมาณค่าแบบจุด

การประมาณค่าแบบจุดสามารถทำได้โดยใช้สถิติพื้นฐานสรุปแนวโน้มสู่ส่วนกลางของค่าพารามิเตอร์ จากการแจกแจงความน่าจะเป็นภายหลัง เช่น ค่าเฉลี่ย มัธยฐาน และฐานนิยม นอกจากนี้ยังสามารถคำนวณค่าส่วนเบี่ยงเบนมาตรฐานเพื่อประเมินระดับความน่าเชื่อถือของค่าประมาณแบบจุดดังกล่าวได้อีกด้วย จากการแจกแจงความน่าจะเป็นภายหลังข้างต้น จะได้ว่า

# mean = expected value of theta
mean.theta<-sum(post.theta*theta)
sd.theta<-sqrt(sum(post.theta*(theta-mean.theta)^2))
mean.theta
## [1] 0.5540441
sd.theta
## [1] 0.1145753

(2) การประมาณค่าแบบช่วง

ช่วงการประมาณที่ได้จากวิธีแบบเบส์จะเรียกว่า ช่วงความน่าเชื่อถือ (credible interval) การประมาณค่าพารามิเตอร์แบบช่วงจากการแจกแจงความน่าจะเป็นภายหลังสามารถทำได้หลายวิธีการ วิธีการง่าย ๆ คือการใช้ค่าเฉลี่ยบวกลบด้วยส่วนเบี่ยงเบนมาตรฐานของพารามิเตอร์ อีกวิธีการหนึ่งที่มีประสิทธิภาพมากกว่าคือการหาช่วงแบบ highest density interval (HDI) กล่าวคือเป็นช่วงการประมาณที่มีความหนาแน่นมากที่สุดที่ทำให้ค่าความน่าจะเป็นของพารามิเตอร์ภายในช่วงดังกล่าวเท่ากับค่าที่กำหนดเช่น ช่วง 95% HDI ของพารามิเตอร์ \(\theta\) ในข้างต้นจะมีค่าประมาณ \([0.4,0.8]\)

ช่วงดังกล่าวมีความแตกต่างจาก 95% confidence interval อย่างไร??

(3) การทดสอบสมมุติฐาน

การทดสอบสมมุติฐานด้วยวิธีการแบบเบส์มีความยืดหยุ่นสูงมาก และสามารถทำได้หลายวิธีการ วิธีการแรกเรียกว่า maximum a posteriori (MAP) test มีรายละเอียดดังนี้

กำหนดให้ \(H_0\) และ \(H_1\) เป็นสมมุติฐานที่ต้องการเปรียบเทียบ ในกรณีนี้จะเห็นว่ามีความไม่แน่นอนว่าสมมุติฐานใดเป็นสมมุติฐานที่ถูกต้องกันแน่ ดังนั้นจึงสามารถกำหนดการแจกแจงความน่าจะเป็นก่อนหน้าของสมมุติฐานทั้งสองได้ ดังนี้

\(p(H_0) = p_0\) และ \(p(H_1)=p_1\) โดยที่ \(p_0+p_1=1\)

และกำหนดให้ \(p(\bf{y}|H_0)\) กับ \(p(\bf{y}|H_1)\) คือฟังก์ชันภาวะความควรจะเป็นบนสมมุติฐานทั้งสอง ซึ่งจากทฤษฎีบทของเบส์จะได้ว่า

\(p(H_0|\bf{y})=\frac{p(\bf{y}|H_0)p(H_0)}{p(\bf{y})}\)

\(p(H_1|\bf{y})=\frac{p(\bf{y}|H_1)p(H_1)}{p(\bf{y})}\)

การตัดสินใจว่าควรเลือกเชื่อสมมุติฐานใด สามารถทำได้โดยดูจากอัตราส่วนที่เรียกว่า posterior odd ซึ่งมีค่าเท่ากับ

\(\frac{p(H_0|\bf{y})}{p(H_1|\bf{y})}=\frac{p(\bf{y}|H_0)p(H_0)}{p(\bf{y}|H_1)p(H_1)}\)

โดยหากอัตราส่วนข้างต้นมีค่ามากกว่า 1 แสดงว่า สมมุติฐาน \(H_0\) มีแนวโน้มน่าเชื่อถือมากกว่า ในทางกลับกันหากอัตราส่วนดังกล่าวมีค่าน้อยกว่า 1 แสดงว่าสมมุติฐาน \(H_1\) มีแนวโน้มน่าเชื่อถือมากกว่า

จากตัวอย่างโยนเหรียญ หากต้องการทดสอบว่าเหรียญมีความเที่ยงตรงหรือไม่ อาจกำหนดสมมุติฐานเป็นดังนี้

\(H_0: \theta = 0.5\) vs \(H_1: \theta>0.5\)

prior.H<-c(0.5,0.5) # prior for H0 and H1 respectively.
# calculate likelihood
L.H0<-0.5^6*(1-0.5)^4
theta.H1<-theta[theta>0.5]
L.H1<-sum(theta.H1^6*(1-theta.H1)^4)

# calculate posterior odd
(prior.H[1]*L.H0)/(prior.H[2]*L.H1)
## [1] 0.3727444

อีกวิธีการหนึ่งที่สามารถทำได้เรียกว่า minimum cost hypothesis test รายละเอียดมีดังนี้

กำหนดให้

  • \(C_{10}\) คือ cost ของการเลือก \(H1\) โดยที่ \(H_0\) ถูกต้อง (cost ของการเกิด type I error)
  • \(C_{01}\) คือ cost ของการเลือก \(H_0\) โดยที่ \(H_1\) ถูกต้อง (cost ของการเกิด type II error)

เกณฑ์การพิจารณาจะใช้ค่าความเสี่ยงภายหลัง (posterior risk) ดังนี้

\(\frac{p(H_0|\bf{y})C_{10}}{p(H_1|\bf{y})C_{01}}\)

จากตัวอย่างโยนเหรียญสมมุติว่า \(C_{10}=5C_{01}\)

c10<-5
c01<-1

# posterior risk of accepting H0
(prior.H[2]*L.H1)*c01
## [1] 0.001309962
# posterior risk of accepting H1
(prior.H[1]*L.H0)*c10
## [1] 0.002441406
# risk ratio
(prior.H[1]*L.H0)*c10/(prior.H[2]*L.H1)*c01
## [1] 1.863722

อีกวิธีการที่ง่ายกว่าสองวิธีการแรกมากคือ วิธี ROPE ซึ่งเป็นการเปรียบเทียบกันระหว่างช่วงความน่าเชื่อถือแบบ HDI กับช่วง ROPE ที่เรียกว่า region of practical equivalence ซึ่งเป็นช่วงที่ยอมรับได้ในทางปฏิบัติ (โดยผู้วิเคราะห์) ว่าหากพารามิเตอร์ยังมีค่าอยู่ในช่วง ROPE ดังกล่าวอยู่ แสดงว่ายังไม่แตกต่างจาก \(H_0\) โดยหากช่วง HDI ครอบคลุมช่วง ROPE เราจะสามารถยอมรับ \(H_0\) ได้ แต่หากช่วง HDI ไม่มีส่วนใดที่ครอบคลุมช่วง ROPE เลย แสดงว่าไม่ช่วงเชื่อถือ \(H_0\) อย่างยิ่ง

การคัดเลือกโมเดล (model selection)

การคัดเลือกโมเดลเป็นวิธีการหนึ่งที่สามารถใช้ในการทดสอบมสมมุติฐานของผู้วิเคราะห์ รวมทั้งใช้เปรียบเทียบและคัดเลือกโมเดลที่เหมาะสมที่สุดสำหรับอธิบายข้อมูล สถิติตัวหนึ่งที่มีความสำคัญคือ bayes factor วัตถุประสงค์ของ bayes factor ถูกใช้เพื่อวัดว่าข้อมูลตัวอย่างที่เก็บรวบรวมมมานี้สนับสนุนโมเดลหรือสมมุติฐานตัวใดมากกว่ากัน bayes factor มีค่าเท่ากับอัตราส่วนระหว่างค่าภาวะความควรจะเป็นของสมมุติฐานสองตัวบนชุดข้อมูลตัวอย่างเดียวกัน สมมุติให้ \(H_1\) คือโมเดลตามสมมุติฐานที่ 1 และ \(H_2\) คือโมเดลตามสมมุติฐานที่ 2 และ \(\bf{x}\) คือชุดข้อมูลตัวอย่าง จะได้ว่าสถิติ bayes factor (BF) มีค่าเท่ากับ

\(BF = \frac{p(\bf{x}|H_1)}{p(\bf{x}|H_2)}=\frac{p(H_1|\bf{x})}{p(H_2|\bf{x})}\times \frac{p(H_2)}{p(H_1)}\)

จากนิยามข้างต้นจะเห็นว่า bayes factor ตามนิยามดังกล่าวสามารถเขียนใหม่ในรูปผลคูณระหว่างอัตราส่วนสองตัว ได้แก่ อัตราส่วนระหว่างการแจกแจงความน่าจะเป็นภายหลังกับอัตราส่วนระหว่างการแจกแจงความน่าจะเป็นก่อนหน้าของสมมุติฐานที่ 1 ต่อสมมุติฐานที่ 2

จากนิยามนี้สามารถแปลผลได้ว่าหาก \(BF\) มีค่ามากกว่า 1 นั่นหมายถึงข้อมูลที่มีสนับสนุนโมเดลตามสมมุติฐาน \(H_1\) มากกว่า ในทางกลับกันหาก \(BF\) มีค่าน้อยกว่า 1 นั่นหมายความว่าข้อมูลที่มีสนับสนุนโมเดลตามสมมุติฐาน \(H_2\) มากกว่า Harold Jeffery (1998) ได้เสนอสเกลการแปลความหมายสถิติ bayes factor ไว้ดังนี้

การคำนวณและการใช้งาน bayes factor จะกล่าวถึงในสัปดาห์ที่สองของหัวข้อนี้

นอกจากสถิติ bayes factor แล้วยังมีสถิติอีกหลายตัวที่สามารถใช้เปรียบเทียบและคัดเลือกโมเดลแบบเบส์ได้ เช่น deviance information criterion (DIC), widely applicable information criterion (WAIC) และ leave-one-out cross-validation (LOO) สถิติทั้งหมดที่กล่าวมานี้ใช้สะท้อนความแตกต่างกันระหว่างข้อมูลตัวอย่าง \(y\) กับค่าทำนายที่ generate จากโมเดลการวิเคราะห์ \(\hat{y}\) ค่าที่ได้จากสถิติดังกล่าวจึงใช้สะท้อนความคลาดเคลื่อนของแต่ละโมเดลเมื่อเปรียบเทียบกับข้อมูลจริง โมเดลที่มีค่าสถิติดังกล่าวต่ำกว่าจะเป็นโมเดลที่มีความเหมาะสมมากกว่า

Monte Carlo Markov Chain (MCMC)

ตัวอย่างที่แสดงให้ดูในหัวข้อข้างต้นเป็นตัวอย่างที่ง่ายมากทำให้การคำนวณการแจกแจงความน่าจะเป็นภายหลังสามารถคำนวณมือได้โดยสะดวก นอกจากนี้ยังสามารถพิสูจน์ในรูปของสูตรปิดทางคณิตศาสตร์ได้ด้วย อย่างไรก็ตามในสถานการณ์ทั่วไปการคำนวณการแจกแจงความน่าจะเป็นภายหลังดังกล่าวมักมีความซับซ้อนและไม่สามารถคำนวณมือได้ โดยเฉพาะในโมเดลพหุระดับ ปัจจุบันวิธีการที่ถูกนำมาแทนการคำนวณทางคณิตศาสตร์คือการใช้ลูกโซ่มาร์คอฟมอนติคาร์โล (MCMC) วิธีการนี้ถูกออกแบบให้สุ่มตัวอย่างพารามิเตอร์ภายใต้กระบวนการลูกโซ่ของมาร์คอฟ ซึ่งรับประกันว่าเมื่อจำนวนรอบของการสุ่มตัวอย่างมีจำนวนมากเพียงพอ ตัวอย่างของพารามิเตอร์ดังกล่าวจะถูกสุ่มมาจากการแจกแจงของประชากรที่เป็นการแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ตามต้องการ ในการทำงานผู้วิเคราะห์จะใช้ตัวอย่างสุ่มที่ถูกตรวจสอบคุณสมบัติแล้ว มาประมาณการแจกแจงความน่าจะเป็นภายหลัง การอนุมานเชิงสถิติต่าง ๆ ที่ได้กล่าวไว้บ้างก่อนหน้านี้สามารถทำได้โดยตรงผ่านการวิเคราะห์ตัวอย่างของพารามิเตอร์ดังกล่าว

ปัจจุบันโปรแกรมสำเร็จรูปที่ช่วยทำ MCMC ให้กับผู้วิเคราะห์มีหลายตัว ได้แก่ BUGS, JAGS, Stan รวมทั้ง R อีกหลาย package นอกจากนี้โปรแกรม Mplus และ MLWins ยังมีโมดูลสำเร็จรูปสำหรับประมาณค่าพารามิเตอร์แบบเบส์สำหรับโมเดลพหุระดับอีกด้วย

บทเรียนนี้จะทำ MCMC ด้วย 2 วิธีการ วิธีการแรกคือการทำ MCMC บนโปรแกรม JAGS ที่สั่งการทำงานบนโปรแกรม R และวิธีการที่สองคือการทำ MCMC ด้วย package-brms ซึ่งเป็น high-level API ของโปรแกรม Stan โดยถูกพัฒนาขึ้นมาสำหรับวิเคราะห์โมเดลพหุระดับโดยเฉพาะ จุดเด่นของ package-brms คือใช้หลักภาษาที่ง่ายเกือบจะเหมือนกับ package-lme4 ที่ได้เรียนไปในบทเรียนก่อนหน้า นอกจากนี้ยังสามารถวิเคราะห์ได้หลายโมเดลตั้งแต่โมเดลแบบตัวแปรตามตัวเดียวไปถึงตัวแปรตามหลายตัว และมีขอบเขตครอบคลุมไปถึง generalized linear model แบบพหุระดับอีกด้วย จุดเด่นอีกประการหนึ่งของ package-brms คืออิงกับโปรแกรม Stan ที่ใช้อัลกอริทึม Hamiltonian Monte Carlo ร่วมกับอัลกอริึม No-U-Turn Sampler (NUTS) ซึ่งมีประสิทธิภาพสูงกว่า Gibb sampler ที่ใช้ใน BUGS และ JAGS อัลกอริทึมดังกล่าวจะลู่เข้าสู่สถานะคงที่ (steady state) ซึ่งก็คือได้ตัวอย่างที่สุ่มจากการแจกแจงความน่าจะเป็นภายหลังไวจึงเหมาะกับโมเดลซับซ้อนที่มีพารามิเตอร์จำนวนมาก นอกจากนี้อัลกอริทึมดังกล่าวยังไม่ต้องการการกำหนดการแจกแจงความน่าจะเป็นก่อนหน้าแบบวงศ์คู่สังยุค (conjugacy prior) เหมือนกับอัลกอริทึม Gibb sampler อีกด้วย จึงทำให้มีความยืดหยุ่นในการวิเคราะห์ที่สูงขึ้นมาก อย่างไรก็ตามข้อจำกัดของ Stan คือการสุ่มตัวอย่างในแต่ละรอบจะใช้เวลาที่มากกว่า BUGS และ JAGS (Bürkner, 2017)

JAGS (Just Another Gibb Sampler)

JAGS สามารถรันได้บน 3 platforms ได้แก่ Mac, Windows และ Linux บทเรียนนี้จะควบคุมโปรแกรม JAGS ผ่านโปรแกรม R โดยใช้ package-rjags ดังนั้นผู้วิเคราะห์จะใช้โปรแกรม R ในการทำงานส่วนใหญ่ ท้ังการเตรียมข้อมูล และการวิเคราะห์ผลลัพธ์จากลูกโซ่มาร์คอฟ ส่วนที่จะเขียนเป็นภาษาของ JAGS คือส่วนการระบุโมเดลการวิเคราะห์เท่านั้น การกำหนดโมเดลใน JAGS ใข้หลักภาษาเดียวกับ BUGS โดยทั่วไปโมเดลการวิเคราะห์แบบเบส์ในโปรแกรม JAGS จะมีส่วนประกอบ 3 ส่วน ได้แก่

  1. likelihood function(s)

  2. prior distribution

  3. deterministic model

prior และ likelihood function เป็นส่วนประกอบที่เรียกว่า stochastic node ซึ่งนิยามโดยใช้สัญลักษณ์ ~ ส่วน deterministic model เป็นส่วนค่าคงที่หรือส่วนโมเดลเชิงคณิตศาสตร์ที่ใช้แสดงความสัมพันธ์ระหว่างตัวแปรการนิยามส่วนนี้จะใช้สัญลักษณ์ <-

ตัวอย่างด้านล่างแสดงการเขียนโมเดลการวิเคราะห์การถดถอยอย่างง่ายแบบเบส์ เพื่อวิเคราะห์ความสัมพันธ์ระหว่างผลสัมฤทธิ์ทางการเรียน (ACH) กับความตั้งใจเรียนของนักเรียน (ATT)

กำหนดให้ \(y\) แทน ACH และ \(x\) แทน ATT เนื่องจากตัวแปรตามในโมเดลเป็นตัวแปรต่อเนื่อง ในกรณีนี้จึงกำหนดให้

\(y_i \sim N(\mu_i, \sigma^2)\)

สังเกตว่าการแจกแจงของค่าสังเกตในข้างต้นมีค่าเฉลี่ยที่ห้อย \(i\) ทั้งนี้เป็นเพราะเรากำลังวิเคราะห์การถดถอย ดังนั้นค่าเฉลี่ยของ \(y\) จึงมีการเปลี่ยนแปลงไปตามค่าของตัวแปรอิสระ \(x_i\) ค่าเฉลี่ยดังกล่าวจึงเป็นค่าเฉลี่ยที่มีเงื่อนไข และสามารถเขียนเป็นสมการแสดงความสัมพันธ์ได้ดังนี้

\(\mu_i = \beta_0 + \beta_1x_i\)

โมเดลข้างต้นมีพารามิเตอร์จำนวน 3 ตัว จำแนกเป็นสองกลุ่ม กลุ่มแรกคือพารามิเตอร์อิทธิพลคงที่ (fixed-effect parameter) ได้แก่ \(\beta_0\) และ \(\beta_1\) และกลุ่มที่สองคือพารามิเตอร์ความแปรปรวนของอิทธิพลสุ่ม ได้แก่ \(\sigma^2\)

การแจกแจงความน่าจะเป็นก่อนหน้า

การแจกแจงความน่าจะเป็นก่อนหน้าใช้สะท้อนความเชื่อเบื้องต้นในพารามิเตอร์ต่าง ๆ ของโมเดล ความเชื่อดังกล่าวมาได้จากทั้งผลการศึกษาในอดีต หลักเหตุผล และประสบการณ์ของนักวิจัย นอกจากนี้ยังควรต้องพิจารณาธรรมชาติหรือค่าที่เป็นไปได้ของพารามิเตอร์ประกอบด้วย

โดยทั่วไปการกำหนดการแจกแจงความน่าจะเป็นก่อนหน้าสามารถทำได้ 2 ลักษณะ ลักษณะแรกเรียกว่า การแจกแจงความน่าจะเป็นก่อนหน้าแบบไม่ให้สารสนเทศ (non-informative prior distribution) และการแจกแจงความน่าจะเป็นก่อนหน้าแบบให้สารสนเทศ (informative prior distribution)

การแจกแจงความน่าจะเป็นก่อนหน้าของพารามิเตอร์อิทธิพลคงที่สามารถเลือกใช้ได้หลายการแจกแจง เช่น การแจกแจงแบบ uniform ที่เป็น noninformative prior การแจกแจงแบบปกติ ที่เป็นไปได้ทั้ง noninformative และ informative prior ขึ้นอยู่กับการกำหนดพารามิเตอร์ค่าเฉลี่ย (mean) และค่าความถูกต้อง (precision) โดยที่ค่าความถูกต้องดังกล่าวเป็นส่วนกลับของพารามิเตอร์ความแปรปรวน (variance) ดังนี้ precision = 1/variance

จากตัวอย่างข้างต้นกำหนดให้

\(\beta_0 \sim N(0,10^{-4})\)

\(\beta_1 \sim N(0,10^{-4})\)

การกำหนดข้างต้นคือการแจกแจงความน่าจะเป็นแบบปกติที่มีค่าเฉลี่ยเท่ากับ 0 และความแปรปรวนเท่ากับ 10000 ซึ่งเรียกได้ว่าเป็นการแจกแจงความน่าจะเป็นก่อนหน้าแบบไม่ให้สารสนเทศ รูปด้่านล่างแสดงลักษณะการแจกแจงแบบปกติข้างต้น จากรูปจะเห็นว่าการแจกแจงดังกล่าวมีขอบเขตที่กว้างพารามิเตอร์ส่วนใหญ่จะมีค่าตกอยู่ในช่วง -200 ถึง 200

การแจกแจงความน่าจะเป็นก่อนหน้าของส่วนเบี่ยงเบนมาตรฐานหรือความแปรปรวน มีความจำเป็นต้องเลือกการแจกแจงที่ domain ของการแจกแจงมีค่าไม่ติดลบ โดยปกติมักกำหนดให้ความแปรปรวนดังกล่าวมีการแจกแจงแบบแกมมาผกผัน (inverse gamma distribution)

จากตัวอย่างข้างต้น กำหนดให้

\(sigma^2 \sim IG(a,b)\)

โดยที่ \(IG\) คือการแจกแจงแบบ inverse-gamma

Note: ในเชิงทฤษฎีหากกำหนดให้ \(sigma^2 \sim IG(a,b)\) จะได้ว่า

\(p(\sigma^2|a,b) \propto (\sigma^2)^{-a-1}exp(\frac{-b}{\sigma^2})\)

หากกำหนดให้ \(a \rightarrow 0\) และ \(b \rightarrow 0\) แล้วฟังก์ชันความหนาแน่นในข้างต้นจะเขียนใหม่ได้เป็น

\(p(\sigma^2|a,b) \propto \frac{1}{\sigma^2}\) เรียกว่า Jeffery prior ซึ่งเป็น noninformative prior ที่ใช้ได้ตัวหนึ่งของพารามิเตอร์ความแปรปรวน

ดังนั้นในกรณีนี้จึงกำหนดให้

\(sigma^2 \sim IG(0.0001,0.0001)\)

รูปด้านล่างแสดงรายการของการแจกแจงความน่าจะเป็นที่สามารถกำหนดได้ใน BUGS และ JAGS

จากการระบุโมเดลข้างต้นทำให้สามารถเขียน syntax เพื่อระบุโมเดลใน JAGS ได้ดังนี้

model{
  
  for(i in 1:n)
  {
    # likelihood function
    y[i]~dnorm(mu[i],tau) #where tau is precision = 1/sigma2
    
    # deteministic model
    mu[i]<-b0+b1x[i] 
  }
  
  # prior distribution
  b0~dnorm(0,10^-4)
  b1~dnorm(0,10^-4)
  
  tau~gamma(0.0001,0.0001)
  sigma2<-1/tau # deteministic model
}

การประมวลผลด้วย JAGS

การประมวลผลด้วย JAGS มี 4 ขั้นตอนหลักได้แก่

(1) เตรียมข้อมูล

(2) ระบุโมเดล

(3) ส่งข้อมูลและโมเดลไปประมวลผลบน JAGS

(4) ตรวจสอบการลู่เข้าของ MCMC

(5) วิเคราะห์และแปลผล

การระบุโมเดล

ขั้นตอนนี้สามารถทำได้ดังตัวอย่างข้างต้น และเมื่อระบุโมเดลในโปรแกรม R เรียบร้อยแล้ว ผู้วิเคราะห์จำเป็นต้องเขียนโมเดลดังกล่าวในรูปของไฟล์ .txt โดยใช้คำสั่ง writeLines(modelString, con="model.txt") ดังตัวอย่างด้านล่าง

modelString<-"model{
  
  for(i in 1:n)
  {
    # likelihood function
    y[i]~dnorm(mu[i],tau) #where tau is precision = 1/sigma2
    
    # deteministic model
    mu[i]<-b0+b1*x[i] 
  }
  
  # prior distribution
  b0~dnorm(0,10^-4)
  b1~dnorm(0,10^-4)
  
  tau~dgamma(0.0001,0.0001)
  sigma2<-1/tau # deteministic model
}" #Wrap model syntax into String object

writeLines(modelString, con="/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/simReg.txt") # write modelString to model.txt

การสั่งประมวลผล

สมมุติว่าข้อมูลที่ใช้ในการวิเคราะห์เป็นดังรูป โดยเก็บไว้ในตัวแปร dat.reg (ข้อมูลในตัวอย่างนี้ได้จากการจำลองด้วยเทคนิค Monte Carlo) โดยใช้ syntax ด้านล่าง

ATT<-rnorm(100,30,10)
ACH<-30+0.55*ATT+rnorm(100,0,5)
dat.reg<-data.frame(ATT,ACH)
head(dat.reg)
##         ATT      ACH
## 1 32.133899 52.34851
## 2  8.208871 37.08404
## 3 37.542008 54.13731
## 4 13.544627 38.93896
## 5 35.487326 44.67005
## 6 42.563916 50.66076
par(mar=c(4,4,1,1))
plot(ATT,ACH, pch=16, xlab="ATT", ylab="ACH")

การประมวลผลโมเดลข้างต้นผ่านโปรแกรม JAGS จะต้องใช้คำสั่ง 2 ตัวได้แก่ฟังก์ชัน jag.model() และ coda.samples() ที่มีอากิวเมนท์สำคัญดังนี้

setwd("/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/")

library(rjags)
dataList<-list(y=dat.reg$ACH, x=dat.reg$ATT, n=length(dat.reg$ACH))
initsList<-list(b0=1, b1=1, tau=1)

model<-jags.model(file = "simReg.txt",
                  data = dataList,
                  inits = initsList,
                  n.chains = 3)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 100
##    Unobserved stochastic nodes: 3
##    Total graph size: 412
## 
## Initializing model
update(model, n.iter=2000) # n.burnin
samples<-coda.samples(model,
                      variable.names=c("b0","b1","sigma2"),
                      n.iter = 10000,
                      thin=1)
head(samples)
## [[1]]
## Markov Chain Monte Carlo (MCMC) output:
## Start = 2001 
## End = 2007 
## Thinning interval = 1 
##            b0        b1   sigma2
## [1,] 30.78025 0.5202950 21.90599
## [2,] 30.60087 0.4840719 24.14094
## [3,] 31.66543 0.4708754 37.18908
## [4,] 31.32383 0.4505433 21.38748
## [5,] 31.89650 0.4600205 24.17811
## [6,] 31.78018 0.4470285 30.98406
## [7,] 32.44455 0.4417329 24.82733
## 
## [[2]]
## Markov Chain Monte Carlo (MCMC) output:
## Start = 2001 
## End = 2007 
## Thinning interval = 1 
##            b0        b1   sigma2
## [1,] 29.98931 0.4762332 30.62127
## [2,] 29.89494 0.5395243 21.85398
## [3,] 29.32648 0.5483189 30.41757
## [4,] 29.36315 0.5429597 31.03104
## [5,] 29.14981 0.5379856 22.40325
## [6,] 30.14993 0.5327550 26.61984
## [7,] 29.80438 0.5197910 23.47122
## 
## [[3]]
## Markov Chain Monte Carlo (MCMC) output:
## Start = 2001 
## End = 2007 
## Thinning interval = 1 
##            b0        b1   sigma2
## [1,] 31.81481 0.4581798 24.36047
## [2,] 31.77167 0.4908677 20.95915
## [3,] 31.31646 0.4657395 25.43489
## [4,] 31.52530 0.4829907 30.82938
## [5,] 31.05918 0.4862660 22.68679
## [6,] 30.97493 0.5190658 24.36354
## [7,] 30.30192 0.5135753 24.93864
## 
## attr(,"class")
## [1] "mcmc.list"

อีก package หนึ่งที่สามารถเรียก JAGS จาก R ได้คือ package-runjags ซึ่งมีจุดเด่นคือสามารถสั่งประมวลแบบคู่ขนาน (parallel) แยกตาม core ของ CPU ได้ในกรณีที่ผู้วิเคราะห์กำหนดให้สุ่มตัวอย่างจากลูกโซ่มาร์คอฟที่มีจำนวนมากกว่า 1 ลูกโซ่ ตัวอย่างคำสั่งเป็นดังนี้

setwd("/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/")

library(runjags)
dataList<-list(y=dat.reg$ACH, x=dat.reg$ATT, n=length(dat.reg$ACH))
initsList<-list(b0=1, b1=1, tau=1)

runjag.output<-run.jags(method="parallel",
                model="simReg.txt",
                monitor=c("b0","b1","sigma2"),
                data=dataList,
                inits=initsList,
                n.chains=3,
                burnin=2000,
                sample=10000, #n.iter
                thin=1,
                summarise=FALSE,
                plots=FALSE)
  
samples<-as.mcmc.list(runjag.output)
head(samples)
plot(samples)

การวินิจฉัยการลู่เข้าของ MCMC

การวิเคราะห์ผลลัพธ์ที่ได้จาก MCMC มีวัตถุประสงค์ 2 ประการ ประการแรกคือการวิเคราะห์เพื่อตรวจสอบ/วินิจฉัยว่าตัวอย่างที่ได้จาก MCMC ในข้างต้น เป็นตัวอย่างที่สุ่มได้จากการแจกแจงความน่าจะเป็นภายหลังที่เป็นเป้าหมายแล้วหรือไม่ (ลูกโซ่ที่สร้างขึ้นลู่เข้าสู่ posterior distribution เป้าหมายแล้วหรือไม่) และประการที่สองคือการวิเคราะห์ผลลัพธ์จาก MCMC เพื่อตอบวัตถุประสงค์ของการวิเคราะห์ข้อมูล

โดยปกติการพิจารณาการลู่เข้าของ MCMC อาจพิจารณาใน 2 มิติ มิติแรกคือการตรวจสอบว่ามีตัวอย่างส่วนใดส่วนหนึ่งของลูกโซ่ที่ลู่ออกหรือมีค่าอยู่นอกเหนือช่วงที่ควรจะเป็นการแจกแจงความน่าจะเป็นภายหลังหรือไม่ โดยปกติตัวอย่างสุ่มที่ได้จาก MCMC จะมีแนวโน้มที่ลู่เข้าหาการแจกแจงความน่าจะเป็นภายหลังมากขึ้นเรื่อย ๆ แต่ในบางสถานการณ์ที่โมเดลมีความซับซ้อน และการระบุโมเดลทำได้ไม่ดีนัก อาจทำให้การลู่เข้าดังกล่าวเป็นไปได้ช้าถึงช้ามาก หรืออาจไม่สามารถลู่เข้าสู่การแจกแจงความน่าจะเป็นภายหลังที่ต้องการได้เลย

ในความเป็นจริงเราไม่สามารถทราบได้ว่าสถานะคงตัวของลูกโซ่มาร์คอฟนั้น เป็นการแจกแจงความน่าจะเป็นภายหลังจริง ๆ หรือไม่ การพิจารณาจึงต้องอาศัยทั้งวิธีการทางสถิติที่ช่วยตรวจสอบการลู่ออกของลูกโซ่ และวิจารณญาณของผู้วิเคราะห์ ในหัวข้อนี้จะกล่าวถึงวิธีการทางสถิติที่สำคัญบางตัว เช่น

(1) Trace plot พล็อตค่าของพารามิเตอร์ที่สุ่มมาจากอัลกอริทึม MCMC เรียงตามลำดับตั้งแต่ลำดับแรกไปจนถึงลำดับสุดท้าย ตัวอย่างดังรูปด้านล่าง

library(MCMCvis)
MCMCtrace(samples, params=c("b0","b1","sigma2"), pdf=F)

(2) Potential Scale Reduction (\(\hat{R}\) หรือ RSRF)

วิธีหนึ่งที่สามารถใช้ตรวจสอบการลู่ออกของลูกโซ่ คือการสร้างลูกโซ่เพื่อประมาณค่าพารามิเตอร์เดียวกันหลาย ๆ ตัว โดยกำหนดค่าเริ่มต้นของแต่ละลูกโซ่ให้แตกต่างกัน หากลูกโซ่ดังกล่าวมีพฤติกรรมที่แตกต่างกันจะสามารถบ่งชี้ได้ว่าตัวอย่างพารามิเตอร์ที่สุ่มจาก MCMC นี้ลู่ออกจากการแจกแจงความน่าจะเป็นภายหลังที่ต้องการ

การวิเคราะห์นี้สามารถพิจารณาได้จากทั้ง trace plot ในข้างต้น และการใช้ค่าสถิติ \(\hat{R}\) ที่มีคำนวณจากอัตราส่วนระหว่างความแปรปรวนรวมต่อความแปรปรวนภายในลูกโซ่ของตัวอย่างพารามิเตอร์

\(\hat{R}=\sqrt{\frac{Var(\theta)}{W}}\)

เมื่อ \(Var(\theta)=(1-\frac{1}{n}W+\frac{1}{n}B)\), \(W=\frac{1}{m}\sum_{j=1}^mS^2_j\),

\(S_j^2=\frac{1}{n-1}\sum_{i=1}^n(\theta_{ij}-\overline{\theta}_j)^2\) และ \(B=\frac{n}{m-n}\sum_{j=1}^m(\overline{\theta}_j-\overline{\theta}_{..})^2\)

Rule of Thumb: \(\hat{R}>1.1\) บ่งชี้ว่าลูกโซ่ของพารามิเตอร์นั้นมีการลู่ออก

MCMCsummary(samples)
##              mean         sd       2.5%        50%      97.5% Rhat n.eff
## b0     30.0629002 1.61239595 26.9107602 30.0698484 33.2015106    1  1633
## b1      0.5213053 0.05160855  0.4217729  0.5211085  0.6222152    1  1607
## sigma2 26.9798622 3.93882192 20.3574546 26.6128567 35.6623178    1 27719

อีกวิธีการคือการพิจารณาเปรียบเทียบค่าคลาดเคลื่อนมาตรฐานระหว่าง Naive SE กับ Time-series SE หากทั้งสองค่านี้มีความแตกต่างกันมากในพารามิเตอร์ตัวใด บ่งชี้ว่าตัวอย่างของพารามิเตอร์ตัวนั้นเกิดอัตสหสัมพันธ์ซึ่งกันและกันสูง

summary(samples)
## 
## Iterations = 2001:12000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean      SD Naive SE Time-series SE
## b0     30.0629 1.61240 0.009309       0.039919
## b1      0.5213 0.05161 0.000298       0.001288
## sigma2 26.9799 3.93882 0.022741       0.023665
## 
## 2. Quantiles for each variable:
## 
##           2.5%     25%     50%    75%   97.5%
## b0     26.9108 28.9716 30.0698 31.166 33.2015
## b1      0.4218  0.4863  0.5211  0.556  0.6222
## sigma2 20.3575 24.1980 26.6129 29.370 35.6623

การแก้ปัญหาในกรณีนี้มีหลายวิธีการซึ่งต้องเลือกทำตามความเหมาะสม

  • กำหนด burnin period เพื่อตัดตัวอย่างส่วนแรกที่ยังไม่ลู่เข้าสู่การแจกแจงความน่าจะเป็นภายหลังออกจากการวิเคราะห์

  • กำหนดจำนวน iteration เพิ่มขึ้นสำหรับ chain ที่ลู่เข้าช้า

  • บางกรณีอาจต้องทำการระบุค่าเริ่มต้นใหม่ เพื่อให้การลู่เข้าทำได้ง่ายขึ้น

  • Reparameterization

มิติที่สองคือการพิจารณาว่าจำนวนรอบของการทวนซ้ำ (iteration) มีความเพียงพอหรือไม่ ปัจจัยหนึ่งที่ควรนำมาพิจารณาร่วมคืออัตสหสัมพันธ์ระหว่างตัวอย่าง ทั้งนี้เป็นเพราะลูกโซ่มาร์คอฟเป็นกระบวนการที่ตัวอย่างแต่ละตัวที่สุ่มขึ้นมาจะมีความสัมพันธ์กับตัวอย่างก่อนหน้า ดังนั้นหาความสัมพันธ์ดังกล่าวมีค่าสูงมากเกินไปจะทำให้ตัวอย่างที่ได้ยังไม่เป็นตัวแทนของการแจกแจงความน่าจะเป็นภายหลังที่ต้องการ

วิธีการหนึ่งที่ใช้พิจารณาว่าตัวอย่างมีความเพียงพอแล้วหรือไม่ คืออาจพิจารณาได้จากตัวชี้วัด effective sample size (ESS) โดยหากำหนดให้ \(T\) คือจำนวน iteration ของ MCMC ที่สุ่มตัวอย่าง \(x_1, x_2, ..., x_T\) ตามลำดับ และ \(n_0\) คือระยะห่างที่น้อยที่สุดที่ทำให้ \(x_t\) กับ \(x_{t+n_0}\) เป็นอิสระซึ่งกันและกัน (หรือมีความสัมพันธ์กันน้อยมากจนตัดทิ้งได้) ดังนั้นตัวอย่างที่เป็นอิสระซึ่งกันและกัน คือตัวอย่างย่อย (sub-sample)

\(x_{n_0}, x_{2n_0},x_{3n_0}...,x_{T}\)

ดังนั้นค่า ESS จะมีค่าเท่ากับ \(T/n_0\) ซึ่งหมายถึงลูกโซ่ขนาด \(T\) ที่สร้างขึ้นมีประสิทธิภาพจริงเทียบเท่าขนาดตัวอย่างเท่าใด

อีกวิธีการหนึ่งที่สามารถใช้ตรวจสอบอัตสหสัมพันธ์ระหว่างตัวอย่างในลูกโซ่มาร์คอฟได้ อีกทั้งยังสามารถใช้ประมาณระยะห่าง \(n_0\) ได้ คือการพิจารณาค่าสัมประสิทธิ์อัตสหสัมพันธ์ (autocorrelation coefficient) และแผนภาพอัตสหสัมพันธ์ (autocorrelation plot)

autocorr.diag(samples)
##                b0         b1      sigma2
## Lag 0  1.00000000 1.00000000 1.000000000
## Lag 1  0.89674469 0.89729301 0.018382972
## Lag 5  0.58881031 0.58970557 0.002997731
## Lag 10 0.34899078 0.35034313 0.004552666
## Lag 50 0.03019034 0.03073331 0.004531337
autocorr.plot(samples[[1]])

เมื่อเกิดสถานการณ์ที่ตัวอย่างในลูกโซ่มีความสัมพันธ์กันอีกสูง การแก้ปัญหาอาจทำได้โดยตัดตัวอย่างลูกโซ่ระหว่าง \(x_{n_0}\) กับ \(x_{2n_0}\) ออก เรียกว่า thining และในบางกรณีอาจจำเป็นต้องเพิ่มจำนวน iteration เพื่อให้มีจำนวนตัวอย่างเพียงพอที่จะเป็นตัวแทนของการแจกแจงความน่าจะเป็นภายหลังด้วย นอกจากนี้การรวมตัวอย่างจากหลายลูกโซ่เข้าด้วยก็ช่วยแก้ปัญหาดังกล่าวได้

การอนุมานเชิงสถิติแบบเบส์

จากผลการวินิจฉัย MCMC ในข้างต้น เราทำการประมวลผลใหม่ดังนี้

setwd("/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/")

#dataList<-list(y=dat.reg$ACH, x=dat.reg$ATT, n=length(dat.reg$ACH))
#initsList<-list(b0=1, b1=1, tau=1)

model<-jags.model(file = "simReg.txt",
                  data = dataList,
                  inits = initsList,
                  n.chains = 3)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 100
##    Unobserved stochastic nodes: 3
##    Total graph size: 412
## 
## Initializing model
update(model, n.iter=2000) # n.burnin
samples<-coda.samples(model,
                      variable.names=c("b0","b1","sigma2"),
                      n.iter = 30000,
                      thin=15)
MCMCsummary(samples)
##              mean         sd       2.5%        50%      97.5% Rhat n.eff
## b0     29.9358012 1.61122865 26.7793446 29.9515443 33.0316408    1  4086
## b1      0.5251646 0.05136625  0.4260559  0.5245155  0.6249333    1  4067
## sigma2 26.9815883 3.93464366 20.3197188 26.6338984 35.8743312    1  6000
MCMCtrace(samples, params=c("b0","b1","sigma2"), pdf=F)

autocorr.diag(samples)
##                   b0           b1       sigma2
## Lag 0    1.000000000 1.0000000000  1.000000000
## Lag 15   0.189732427 0.1917814786 -0.003137443
## Lag 75   0.026335394 0.0191889064  0.007447446
## Lag 150  0.006470748 0.0103972916 -0.007191825
## Lag 750 -0.003792101 0.0003762561  0.025105044
autocorr.plot(samples[[1]])

การประมาณค่าพารามิเตอร์

summary(samples)
## 
## Iterations = 2015:32000
## Thinning interval = 15 
## Number of chains = 3 
## Sample size per chain = 2000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean      SD  Naive SE Time-series SE
## b0     29.9358 1.61123 0.0208009      0.0252067
## b1      0.5252 0.05137 0.0006631      0.0008054
## sigma2 26.9816 3.93464 0.0507960      0.0507991
## 
## 2. Quantiles for each variable:
## 
##           2.5%   25%     50%    75%   97.5%
## b0     26.7793 28.84 29.9515 31.033 33.0316
## b1      0.4261  0.49  0.5245  0.561  0.6249
## sigma2 20.3197 24.21 26.6339 29.332 35.8743

การประมาณค่าแบบช่วงและการทดสอบสมมุติฐานของพารามิเตอร์

ฟังก์ชัน plotPost() สามารถใช้ประมาณช่วงความน่าเชื่อถือแบบ HDI ได้โดยง่าย ฟังก์ชันนี้มีอาร์กิวเมนท์ได้แก่

  • codasamples

  • cenTend

  • compVal

  • ROPE

  • credMass

\(H_0: \beta_1 =0\)

#library(googledrive)
setwd("/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel")
#drive_download("https://drive.google.com/file/d/0B3lh4V2Mrl14SF94US1YTkNVaHM/view?resourcekey=0-bYI7GCSuypi3a5aFt0O-qg")
source("DBDA2E-utilities.R") # (Kruschke, 2015)
## 
## *********************************************************************
## Kruschke, J. K. (2015). Doing Bayesian Data Analysis, Second Edition:
## A Tutorial with R, JAGS, and Stan. Academic Press / Elsevier.
## *********************************************************************
plotPost(samples[,"b1"], compVal=0, cenTend="median", ROPE=c(0,0.1), credMass=0.95)

##                  ESS      mean    median      mode hdiMass    hdiLow   hdiHigh
## Param. Val. 4069.857 0.5251646 0.5245155 0.5166871    0.95 0.4260469 0.6249327
##             compVal pGtCompVal ROPElow ROPEhigh pLtROPE pInROPE pGtROPE
## Param. Val.       0          1       0      0.1       0       0       1

Note: ROPE ย่อมาจาก region of pratical equivalence ซึ่งก็คือช่วงการยอมรับเชิงปฏิบัติ หลักการง่าย ๆ คือ การทดสอบสมมุติฐาน \(H_0: \theta = \theta_0\) แทบจะมีมีโอกาสที่ค่าเฉลี่ยตัวอย่างของพารามิเตอร์ \(\theta\) ที่สุ่มจาก MCMC จะมีค่าเท่ากับ \(\theta_0\) อย่างพอดี การกำหนดช่วงการยอมรับเชิงปฏิบัติ ทำให้ผู้วิเคราะห์มีเกณฑ์การพิจารณาว่าควรตัดสินใจเชื่อ \(H_0\) หรือไม่ จากการเปรียบเทียบ ROPE กับช่วง HDI หากช่วง HDI ครอบคลุมช่วง ROPE อยู่นั้นแสดงว่ายังมีโอกาสสูงที่พารามิเตอร์จะมีค่าเป็นไปตามสมมุติฐาน \(H_0\)

\(H_0: \beta_1 =0.6\)

plotPost(samples[,"b1"], compVal=0.6, cenTend="median", ROPE=c(0.55,0.65), credMass=0.95)

##                  ESS      mean    median      mode hdiMass    hdiLow   hdiHigh
## Param. Val. 4069.857 0.5251646 0.5245155 0.5166871    0.95 0.4260469 0.6249327
##             compVal pGtCompVal ROPElow ROPEhigh   pLtROPE   pInROPE     pGtROPE
## Param. Val.     0.6 0.07133333    0.55     0.65 0.6831667 0.3091667 0.007666667

ตัวอย่างโยนเหรียญ (รอบที่ 2)

  • จากตัวอย่างที่ 1 จะได้ว่าโมเดลของค่าสังเกตเป็นฟังก์ชันความน่าจะเป็นแบบ Bernuolli ที่มีพารามิเตอร์เป็น \(\theta\) เขียนแทนด้วย
\(y_i \sim Ber(\theta)\) โดยที่ \(i=1,2,...,n\)
  • อย่างไรก็ตามในกรณีนี้ผู้วิเคราะห์ต้องการกำหนดการแจกแจงความน่าจะเป็นก่อนหน้าแบบต่อเนื่องให้กับพารามิเตอร์ความลำเอียง \(\theta\) เนื่องจากพารามิเตอร์ดังกล่าวต้องมีค่าอยู่บนช่วง \([0,1]\) การแจกแจงความน่าจะเป็นที่เหมาะสมจึงเป็นการแจกแจงแบบบีต้า (beta distribution) ที่มีพารามิเตอร์ตำแหน่ง \(a\) และพารามิเตอร์สเกล \(b\) ดังนั้น
\(\theta \sim Beta(a,b)\)

พารามิเตอร์ทั้งสองตัวของการแจกแจงแบบบีต้าทำให้รูปทรงการแจกแจงมีความแตกต่างกัน ซึ่งสามารถใช้สะท้อนความเชื่อเกี่ยวกับความลำเอียงของเหรียญที่แตกต่างกันได้เป็นอย่างดี รูปด้านล่างแสดงตัวอย่างโค้งความหนาแน่นของการแจกแจงแบบบีต้า เมื่อกำหนดพารามิเตอร์ต่าง ๆ

Kruschke (2012)


จากการกำหนดในข้างต้นจะสามารถเขียน syntax ของโมเดลในโปรแรกม JAGS ได้ดังนี้

model{

# likelihood part    
for (i in 1:n)
{
  y[i]~dbern(theta)
  
}
  
# prior part
  theta~beta(a,b)

# deterministic part
  a<-3
  b<-2
}

การประมวลผลโมเดลข้างต้นด้วยโปรแกรม JAGS บน R จำเป็นต้องมีการ wrap up คำสั่งข้างต้นให้เป็น text file เพื่อส่งออกไปยังโปรแกรม R โดยเขียนคำสั่งดังนี้

coin_model="

model{

# likelihood part    
for (i in 1:n)
{
  y[i]~dbern(theta)
  
}
  
# prior part
  theta~dbeta(a,b)

# deterministic part
  a<-3
  b<-2
}

"
writeLines(coin_model, con="/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/coin_model.txt")

ขั้นตอนถัดมาคือการประมาณการแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ \(\theta\) โดยใช้เทคนิค MCMC ผ่านโปรแกรม JAGS ซึ่งสามารถสั่งงานได้ผ่านฟังก์ชัน jag.model() โดยฟังก์ชันนี้มีอาร์กิวเมนท์ที่สำคัญที่เกี่ยวข้องกับการระบุคุณลักษณะของลูกโซ่มาร์คอฟที่ต้องการสร้างขึ้น ได้แก่

  • ข้อมูลค่าสังเกต ในที่นี้สมมุติว่าทำการทดลอง 20 ครั้งได้ผลเป็นดังนี้
y<-rbinom(20,1,0.65)
data.list<-list(y=y, n=20)
  • n.chains
library(rjags)
fit<-jags.model(file="/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/coin_model.txt", data=data.list,
                n.chains = 3)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 20
##    Unobserved stochastic nodes: 1
##    Total graph size: 24
## 
## Initializing model
samples<-coda.samples(fit,variable.names=c("theta"),n.iter=5000)
plot(samples)

autocorr.diag(samples)
##              theta
## Lag 0   1.00000000
## Lag 1   0.01582687
## Lag 5  -0.01269683
## Lag 10  0.01048772
## Lag 50 -0.01040355
summary(samples)
## 
## Iterations = 1:5000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean             SD       Naive SE Time-series SE 
##      0.6796919      0.0908805      0.0007420      0.0007499 
## 
## 2. Quantiles for each variable:
## 
##   2.5%    25%    50%    75%  97.5% 
## 0.4931 0.6189 0.6847 0.7465 0.8412

เหรียญเที่ยงตรงหรือไม่?? —> กำหนด ROPE = [0.48,0.52]

plotPost(samples[,"theta"], compVal=0.5, cenTend="median", ROPE=c(0.48,0.1), credMass=0.95)

##                  ESS      mean    median      mode hdiMass    hdiLow   hdiHigh
## Param. Val. 14532.27 0.6796919 0.6847037 0.7075729    0.95 0.5029832 0.8492614
##             compVal pGtCompVal ROPElow ROPEhigh pLtROPE pInROPE pGtROPE
## Param. Val.     0.5     0.9702    0.48      0.1  0.0186       0       1

Multilevel Modelling

หัวข้อนี้จะกล่าวถึงการวิเคราะห์ multilevel model ด้วยวิธีการแบบเบส์ จากชุดข้อมูล hsb สมมุติว่านักวิจัยมีคำถามวิจัยดังนี้

  1. ผลสัมฤทธิ์ทางการเรียนวิชาคณิตศาสตร์ของนักเรียนมีความแตกต่างกันระหว่างโรงเรียนหรือไม่?

  2. โรงเรียนโดยเฉลี่ยแล้วนักเรียนมีเศราฐานะทางบ้านสูงจะมีค่าเฉลี่ยผลสัมฤทธิ์ทางคณิตศาสตร์สูงด้วยหรือไม่?

  3. ความสัมพันธ์ในระดับนักเรียนระหว่าง SES กับ mathach มีมากน้อยแค่ไหน? ความสัมพันธ์ดังกล่าวมีความแตกต่างกันไปตามโรงเรียนหรือไม่อย่างไร?

  4. เมื่อเปรียบเทียบกันระหว่างโรงเรียนรัฐบาลกับโรงเรียนในเครือคาทอลิก ในมิติของผลสัมฤทธิ์ทางการเรียนวิชาคณิตศาตร์ และระดับความสััมพันธ์ระหว่าง ses กับ mathach มีความแตกต่างกันหรือไม่ เมื่อมีการควบคุมความผันแปรของ mathach ด้วยค่าเฉลี่ยเศรษฐานะของนักเรียนในแต่ละโรงเรียนแล้ว

จงออกแบบการวิเคราะห์เพื่อตอบคำถามวิจัยในข้างต้น

One-Way ANOVA with random effect

null model หรือโมเดลการวิเคราะห์ความแปรปรวนแบบอิทธิพลสุ่มมีสมการดังนี้

level-1 model: \(y_{ij}=\beta_{0j}+\epsilon_{ij}\)

level-2 model: \(\beta_{0j}=\gamma_{00}+u_{0j}\)

combined model: \(y_{ij}=\gamma_{00}+u_{0j}+\epsilon_{ij}\)

โดยที่ตัวแปรตามคือ mathach ในชุดข้อมูล hsb หากดำเนินการวิเคราะห์ข้อมูลด้วยวิธีการแบบดั้งเดิม จะได้ผลการวิเคราะห์เป็นดังนี้

library(foreign)
dat1<-read.spss("/Users/siwachoat/Downloads/hsb1.sav", to.data.frame = TRUE)
dat2<-read.spss("/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/เอกสารประกอบการสอน/2756714/HLM/hsb2.sav", to.data.frame=TRUE)
dat<-merge(dat1,dat2, by="schoolid")
library(lme4)
## Loading required package: Matrix
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
fit<-lmer(mathach~1+(1|schoolid), data=dat1)
ranova(fit) #reduce random-effect term
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## mathach ~ (1 | schoolid)
##                npar logLik   AIC    LRT Df Pr(>Chisq)    
## <none>            3 -23558 47123                         
## (1 | schoolid)    2 -24052 48107 986.12  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(fit)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mathach ~ 1 + (1 | schoolid)
##    Data: dat1
## 
## REML criterion at convergence: 47116.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0631 -0.7539  0.0267  0.7606  2.7426 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  schoolid (Intercept)  8.614   2.935   
##  Residual             39.148   6.257   
## Number of obs: 7185, groups:  schoolid, 160
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  12.6370     0.2444 156.6473   51.71   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
performance::icc(fit)
## # Intraclass Correlation Coefficient
## 
##      Adjusted ICC: 0.180
##   Conditional ICC: 0.180
# reliability of the sample mean (Raudenbush & Bryk (2002) p.46)
nj<-table(dat1$schoolid)
reliaj<-8.614/(8.614+39.148/nj)
hist(reliaj)

mean(reliaj)
## [1] 0.9013778
#calculate level-1 residual
res1<-dat1$mathach-predict(fit)
res1<-data.frame(schoolid=dat1$schoolid, res1)
#calculate level-2 residual
res2<-coef(fit)$schoolid-12.6370
names(res2)<-"res2"
res2<-data.frame(schoolid=as.numeric(names(table(dat1$schoolid))), res2)


#normality
ggplot(res1)+geom_histogram(aes(x=res1))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(res2)+geom_histogram(aes(x=res2))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

library(ggpubr)
ggqqplot(res1$res1)

ggqqplot(res2$res2)

#homogeneity of variance
res1%>%group_by(schoolid)%>%summarise(var.e=var(res1))%>%summary()
##     schoolid        var.e      
##  Min.   :1224   Min.   :12.54  
##  1st Qu.:3018   1st Qu.:31.27  
##  Median :5512   Median :39.94  
##  Mean   :5310   Mean   :39.75  
##  3rd Qu.:7432   3rd Qu.:47.10  
##  Max.   :9586   Max.   :71.93
res1%>%group_by(schoolid)%>%summarise(var.e=var(res1))%>%ggplot(aes(x=var.e))+geom_density()

หากดำเนินการวิเคราะห์ด้วยวิธีการแบบ bayes โดยใช้โปรแกรม JAGS อาจดำเนินการได้ดังนี้

library(rjags)

#generate model file
null.model<-"model{

for (i in 1:n)
{
y[i]~dnorm(mu[i],inv.sigma2)
mu[i]<-beta0j[schoolid[i]]
}

for (j in tab)
{
beta0j[j]~dnorm(gamma00, inv.tau00)
}

inv.sigma2~dgamma(0.0001,0.0001)
inv.tau00~dgamma(0.0001,0.0001)
gamma00~dnorm(0,10^-4) T(0,) #half-normal distribution

sigma2<-inv.sigma2^(-1)
tau00<-inv.tau00^(-1)

icc<-tau00/(tau00+sigma2)

}"

writeLines(null.model,con="null_model.txt")

datlist<-list(y=dat1$mathach,schoolid=dat1$schoolid, n=length(dat1$mathach),
              tab=as.numeric(names(table(dat1$schoolid))))
fit.null<-jags.model(file="null_model.txt",
                     data=datlist,
                     n.chains=3)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 7185
##    Unobserved stochastic nodes: 163
##    Total graph size: 14706
## 
## Initializing model
update(fit.null,n.iter=2000) #n.burnin
output<-coda.samples(fit.null,
                     variable.names=c("gamma00","tau00","sigma2","icc"),
                     n.iter=30000,
                     thin=2)

library(MCMCvis)
MCMCtrace(output, params=c("gamma00","icc"), pdf=F)

MCMCsummary(output)
##               mean         sd       2.5%        50%      97.5% Rhat n.eff
## gamma00 12.6367876 0.24487553 12.1583073 12.6381989 13.1172507    1 45000
## icc      0.1817391 0.01885253  0.1474386  0.1808955  0.2212261    1 40991
## sigma2  39.1627730 0.65410399 37.9016391 39.1535162 40.4598172    1 44151
## tau00    8.7211157 1.10051603  6.7886400  8.6454339 11.0962064    1 40938
summary(output)
## 
## Iterations = 2002:32000
## Thinning interval = 2 
## Number of chains = 3 
## Sample size per chain = 15000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##            Mean      SD  Naive SE Time-series SE
## gamma00 12.6368 0.24488 1.154e-03      1.154e-03
## icc      0.1817 0.01885 8.887e-05      9.312e-05
## sigma2  39.1628 0.65410 3.083e-03      3.114e-03
## tau00    8.7211 1.10052 5.188e-03      5.440e-03
## 
## 2. Quantiles for each variable:
## 
##            2.5%     25%     50%     75%   97.5%
## gamma00 12.1583 12.4714 12.6382 12.8033 13.1173
## icc      0.1474  0.1686  0.1809  0.1939  0.2212
## sigma2  37.9016 38.7157 39.1535 39.6046 40.4598
## tau00    6.7886  7.9512  8.6454  9.4048 11.0962
plotPost(output[,"icc"], compVal=0.05, cenTend="median", ROPE=c(0.04,0.06), credMass=0.95,
         xlab="ICC")

##          ESS      mean    median      mode hdiMass    hdiLow   hdiHigh compVal
## ICC 41391.37 0.1817391 0.1808955 0.1794032    0.95 0.1457078 0.2191428    0.05
##     pGtCompVal ROPElow ROPEhigh pLtROPE pInROPE pGtROPE
## ICC          1    0.04     0.06       0       0       1
library(rjags)

#generate model file
null.model<-"model{

for (i in 1:n)
{
y[i]~dnorm(mu[i],inv.sigma2[schoolid[i]])
mu[i]<-beta0j[schoolid[i]]
}

for (j in tab)
{
beta0j[j]~dnorm(gamma00, inv.tau00)
inv.sigma2[j]~dgamma(0.0001,0.0001)
sigma2[j]<-inv.sigma2[j]^(-1)

}

inv.tau00~dgamma(0.0001,0.0001)
gamma00~dnorm(0,10^-4) T(0,) #half-normal distribution
tau00<-inv.tau00^(-1)
}"

writeLines(null.model,con="null_model.txt")

datlist<-list(y=dat1$mathach,schoolid=dat1$schoolid, n=length(dat1$mathach),
              tab=as.numeric(names(table(dat1$schoolid))))
fit.hetero<-jags.model(file="null_model.txt",
                     data=datlist,
                     n.chains=3)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 7185
##    Unobserved stochastic nodes: 322
##    Total graph size: 15022
## 
## Initializing model
update(fit.hetero,n.iter=2000) #n.burnin
output<-coda.samples(fit.hetero,
                     variable.names=c("gamma00","tau00","sigma2"),
                     n.iter=30000,
                     thin=2)

library(MCMCvis)
MCMCsummary(output)
##                  mean         sd      2.5%       50%     97.5% Rhat n.eff
## gamma00      12.65549  0.2514842 12.160433 12.655888  13.14829    1 44242
## sigma2[1224] 60.19157 13.0677194 39.872462 58.449081  90.65491    1 43915
## sigma2[1288] 53.16940 16.5107789 30.066191 50.173312  93.34146    1 48202
## sigma2[1296] 29.98479  6.5156131 19.869552 29.097814  45.04621    1 45000
## sigma2[1308] 41.80058 15.0983901 21.594343 38.774190  79.14488    1 43913
## sigma2[1317] 31.12595  6.7245490 20.631624 30.210267  46.78136    1 45139
## sigma2[1358] 36.93542 10.4336943 21.831424 35.186547  62.27490    1 44588
## sigma2[1374] 75.43277 22.2474597 43.540262 71.621535 129.38767    1 44960
## sigma2[1433] 16.04133  4.1573219  9.836660 15.405604  25.92200    1 45000
## sigma2[1436] 21.78928  4.9608527 14.178401 21.087452  33.52247    1 45000
## sigma2[1461] 51.46226 13.7561788 31.219932 49.250434  84.58928    1 45000
## sigma2[1462] 41.23397  8.1201423 28.301439 40.212561  59.87144    1 45751
## sigma2[1477] 52.86930  9.8575519 37.028408 51.740710  75.52237    1 45000
## sigma2[1499] 41.90864  8.5942796 28.323192 40.814185  61.53317    1 45000
## sigma2[1637] 55.70535 17.0675853 31.773920 52.783438  97.24919    1 43574
## sigma2[1906] 44.22319  9.0872565 29.905604 43.024882  65.20067    1 45000
## sigma2[1909] 40.72442 12.0302599 23.501970 38.643688  69.78552    1 45689
## sigma2[1942] 32.56230  9.4782120 19.002048 30.980386  55.38963    1 45000
## sigma2[1946] 51.27276 12.3699894 32.584266 49.428249  80.57575    1 45000
## sigma2[2030] 41.28092  8.9668134 27.366108 40.057587  62.26450    1 44440
## sigma2[2208] 38.82988  7.4463692 26.997708 37.934845  55.96154    1 43710
## sigma2[2277] 32.03703  6.0690991 22.234330 31.296887  46.02290    1 45000
## sigma2[2305] 25.89295  4.6239434 18.419459 25.367988  36.54755    1 44896
## sigma2[2336] 35.78740  7.7935331 23.687094 34.738324  53.81243    1 44321
## sigma2[2458] 35.37772  6.9451626 24.317726 34.526280  51.29790    1 45465
## sigma2[2467] 47.85396  9.8813976 32.218197 46.618612  70.90522    1 43913
## sigma2[2526] 23.69674  4.6673615 16.306127 23.134159  34.57381    1 44077
## sigma2[2626] 41.11736 10.1278780 25.901219 39.619390  65.11114    1 44599
## sigma2[2629] 27.64507  5.4027537 19.065542 26.954774  39.99894    1 44623
## sigma2[2639] 36.19415  8.4455039 23.273651 35.045227  56.18398    1 45000
## sigma2[2651] 51.86463 12.6833077 32.668867 49.937409  82.07508    1 44651
## sigma2[2655] 41.53026  8.5117493 28.071720 40.401609  61.29319    1 44513
## sigma2[2658] 33.26988  7.4519265 21.808922 32.246284  50.66791    1 44536
## sigma2[2755] 38.57854  8.3794385 25.506217 37.439966  58.06386    1 45493
## sigma2[2768] 57.43159 18.0613354 32.075967 54.181315 102.02909    1 44939
## sigma2[2771] 47.96107  9.5784090 32.878142 46.737742  70.15289    1 46234
## sigma2[2818] 31.87086  7.3520616 20.544201 30.856273  49.23508    1 43681
## sigma2[2917] 52.78165 12.1475424 34.119049 51.037604  81.29278    1 45385
## sigma2[2990] 19.72811  4.2542055 13.135323 19.173486  29.57098    1 45215
## sigma2[2995] 55.84869 12.3698831 36.606600 54.212405  84.90843    1 44571
## sigma2[3013] 50.60847 10.2515669 34.335236 49.336139  74.20117    1 44502
## sigma2[3020] 41.70495  8.0386673 28.873111 40.713629  60.06811    1 45000
## sigma2[3039] 32.19949 11.4030505 16.956329 29.984724  60.32980    1 45247
## sigma2[3088] 37.93164  9.2291605 23.977790 36.545053  59.81383    1 44641
## sigma2[3152] 56.88409 11.6967789 38.496180 55.370426  83.84153    1 45000
## sigma2[3332] 44.77404 10.9695837 28.151362 43.143046  70.76379    1 44435
## sigma2[3351] 52.18636 12.7072259 32.833498 50.331278  82.03890    1 44321
## sigma2[3377] 51.64357 11.5700697 33.674183 50.079745  78.78349    1 44328
## sigma2[3427] 13.15345  2.8222248  8.736623 12.785176  19.80100    1 45390
## sigma2[3498] 25.97002  5.2919369 17.564168 25.296574  38.18994    1 44511
## sigma2[3499] 31.45782  7.6953345 19.825063 30.341658  49.57740    1 44924
## sigma2[3533] 54.56205 11.7384429 36.166382 53.048400  81.93294    1 45000
## sigma2[3610] 35.87454  6.6020964 25.225146 35.109942  50.95708    1 44821
## sigma2[3657] 62.20887 12.8843325 41.801134 60.576876  91.98888    1 45000
## sigma2[3688] 47.29698 10.8148645 30.614391 45.722289  72.96350    1 45000
## sigma2[3705] 40.55032  9.0918732 26.440599 39.296417  61.97958    1 44033
## sigma2[3716] 75.64338 17.8456689 48.458948 73.001054 117.48385    1 44613
## sigma2[3838] 27.07179  5.4739627 18.448212 26.343909  39.78723    1 45517
## sigma2[3881] 45.88496 10.7714342 29.443579 44.326782  71.49261    1 44448
## sigma2[3967] 49.32437 10.1505469 33.285648 48.019856  72.86548    1 45000
## sigma2[3992] 32.48634  6.6253652 22.033897 31.639625  47.82663    1 45000
## sigma2[3999] 63.68147 14.0315848 41.771119 61.796220  96.44179    1 44880
## sigma2[4042] 42.88707  7.8937138 30.045636 41.997704  60.91079    1 44197
## sigma2[4173] 51.35420 11.5316829 33.611976 49.730882  78.11309    1 45000
## sigma2[4223] 43.00422  9.5897565 28.082518 41.727276  65.60563    1 43404
## sigma2[4253] 33.19243  6.4509104 22.836672 32.397621  47.98253    1 45000
## sigma2[4292] 39.93381  7.2753959 28.219372 39.128877  56.52703    1 45000
## sigma2[4325] 37.66066  7.7088590 25.559442 36.704958  55.60960    1 45000
## sigma2[4350] 55.17794 14.7413463 33.332433 52.847586  90.26246    1 45000
## sigma2[4383] 59.97786 18.6782108 33.838398 56.594689 105.22696    1 44550
## sigma2[4410] 40.21179  9.4435115 25.800708 38.877802  62.31101    1 45000
## sigma2[4420] 45.92915 12.4932123 27.541042 43.901113  75.81138    1 45889
## sigma2[4458] 20.98883  4.5233121 13.885132 20.393734  31.58704    1 45000
## sigma2[4511] 36.62661  7.1252640 25.238617 35.748666  53.01316    1 45359
## sigma2[4523] 32.11103  7.0526311 21.151461 31.177179  48.55198    1 45000
## sigma2[4530] 32.08005  5.9805556 22.532385 31.370734  45.65906    1 45000
## sigma2[4642] 38.35838  7.2570457 26.698959 37.467132  54.92376    1 45819
## sigma2[4868] 31.26781  8.1631762 19.056097 30.049310  50.91516    1 45353
## sigma2[4931] 29.39309  5.6848249 20.333775 28.712575  42.46148    1 45000
## sigma2[5192] 53.85589 15.8712952 31.197746 51.206495  92.15191    1 46876
## sigma2[5404] 37.61934  7.3585860 25.830609 36.738920  54.61382    1 45374
## sigma2[5619] 54.61691  9.9104725 38.674234 53.422470  77.17139    1 44561
## sigma2[5640] 52.25551 10.2331550 35.910943 50.940000  76.04966    1 45000
## sigma2[5650] 43.82994  9.7287861 28.771032 42.482415  66.36990    1 45000
## sigma2[5667] 44.81248  8.4603252 31.335715 43.782942  63.99972    1 45000
## sigma2[5720] 33.67890  6.8762154 22.873839 32.797800  49.50510    1 45000
## sigma2[5761] 44.46468  9.1236501 30.022948 43.347238  65.35517    1 45000
## sigma2[5762] 26.76025  6.7917153 16.587019 25.716064  43.04878    1 45000
## sigma2[5783] 45.67161 13.1608142 26.716212 43.507649  77.71426    1 44609
## sigma2[5815] 34.47724 10.9922264 19.261471 32.499961  61.56169    1 44300
## sigma2[5819] 50.57602 10.6391394 33.959650 49.133487  75.24668    1 45000
## sigma2[5838] 33.30134  9.1853551 19.842431 31.875482  55.33488    1 44572
## sigma2[5937] 26.60941  7.7254317 15.450716 25.331708  45.16086    1 44483
## sigma2[6074] 41.08376  8.1215975 28.156115 40.061159  59.88008    1 44503
## sigma2[6089] 45.04725 12.0059738 27.268417 43.228287  73.90989    1 44506
## sigma2[6144] 44.61165 10.2848589 28.876943 43.153325  68.94483    1 44958
## sigma2[6170] 76.39065 26.8142096 40.705098 71.070778 142.82118    1 45000
## sigma2[6291] 46.19685 11.9690269 28.443755 44.267533  74.81003    1 44036
## sigma2[6366] 41.06237  7.9854519 28.318932 40.062493  59.29037    1 46221
## sigma2[6397] 49.94028  9.4773093 34.827726 48.838631  71.70610    1 44977
## sigma2[6415] 47.60867  9.6312035 32.341571 46.377283  69.74527    1 44451
## sigma2[6443] 37.55948 10.5376479 22.266250 35.815627  62.88573    1 45000
## sigma2[6464] 31.15446  9.0638023 18.119497 29.583505  53.03124    1 45000
## sigma2[6469] 23.67562  4.6462931 16.264451 23.107028  34.35422    1 45000
## sigma2[6484] 58.25340 14.9735380 35.958298 56.000296  93.57001    1 45000
## sigma2[6578] 41.36688  8.1544071 28.339339 40.377228  60.09285    1 46352
## sigma2[6600] 48.59893  9.5948826 33.339704 47.421753  70.60684    1 45000
## sigma2[6808] 51.30961 11.5608818 33.441305 49.707366  78.29894    1 45385
## sigma2[6816] 22.01430  4.3771878 14.997027 21.483280  32.07500    1 45000
## sigma2[6897] 46.08964  9.8717144 30.675178 44.748584  69.23992    1 44560
## sigma2[6990] 29.45524  6.0321837 19.866805 28.686394  43.38169    1 45684
## sigma2[7011] 39.47482 10.5175163 23.875128 37.820573  64.68349    1 45000
## sigma2[7101] 41.85311 12.1840750 24.237028 39.730191  71.20206    1 45472
## sigma2[7172] 33.16096  7.5140550 21.569403 32.120083  50.73435    1 43921
## sigma2[7232] 46.88468  9.6019194 31.788116 45.601123  69.07229    1 44898
## sigma2[7276] 47.24516  9.6098757 32.077532 46.010865  69.51585    1 45406
## sigma2[7332] 33.46157  7.2179354 22.179040 32.487669  50.35167    1 46957
## sigma2[7341] 30.38838  6.3440396 20.468837 29.550022  45.21840    1 44577
## sigma2[7342] 42.33375  8.2046688 29.193700 41.343308  61.14980    1 44135
## sigma2[7345] 54.31211 10.7049913 37.212751 52.963423  79.02026    1 45000
## sigma2[7364] 35.57372  8.0663265 23.198048 34.423508  54.52870    1 44218
## sigma2[7635] 37.92973  7.9663986 25.450623 36.860914  56.38467    1 45202
## sigma2[7688] 21.06221  4.2749509 14.322169 20.519848  30.99096    1 45195
## sigma2[7697] 46.80730 12.6880148 28.146193 44.790827  77.15302    1 45000
## sigma2[7734] 78.51588 26.6764906 42.217812 73.596750 143.72939    1 45969
## sigma2[7890] 40.79607  8.5007120 27.355123 39.711317  60.32602    1 45000
## sigma2[7919] 48.82882 12.1798757 30.532472 46.925516  77.81224    1 44069
## sigma2[8009] 33.34019  7.1742189 22.126490 32.387149  49.96200    1 45490
## sigma2[8150] 26.70163  6.0602000 17.381521 25.843245  40.94066    1 44935
## sigma2[8165] 29.60439  6.3325262 19.665576 28.760315  44.45246    1 45925
## sigma2[8175] 39.59333 10.5725710 24.037960 37.854049  64.92911    1 45362
## sigma2[8188] 47.08454 13.2232550 27.885049 44.850434  79.22685    1 44573
## sigma2[8193] 19.68501  4.5142595 12.728704 19.026518  30.31206    1 44382
## sigma2[8202] 61.96347 15.9745772 38.233317 59.450685 100.16733    1 45825
## sigma2[8357] 48.67640 14.5584921 27.949632 46.159250  84.37082    1 43902
## sigma2[8367] 25.34600 13.0806371 10.705211 22.206306  59.02414    1 42226
## sigma2[8477] 48.94848 12.1509508 30.699126 47.132693  77.74330    1 45000
## sigma2[8531] 58.47321 13.7231198 37.521123 56.591515  90.87648    1 44572
## sigma2[8627] 43.79554  8.8865157 29.665584 42.697393  64.24942    1 45421
## sigma2[8628] 35.29145  6.6543725 24.565609 34.463447  50.50385    1 45000
## sigma2[8707] 43.13820  9.2183214 28.680666 41.952435  64.43249    1 45000
## sigma2[8775] 33.19618  7.1760534 21.911244 32.285265  49.80626    1 46404
## sigma2[8800] 44.87803 12.2748803 26.928511 42.876265  74.53022    1 46530
## sigma2[8854] 31.93908  8.8895217 18.965028 30.478818  53.38356    1 45000
## sigma2[8857] 29.73065  5.4644029 20.897571 29.105310  42.34175    1 45000
## sigma2[8874] 50.67831 12.8735408 31.455818 48.650567  81.40370    1 45122
## sigma2[8946] 43.97581  8.5301792 30.355238 42.958567  63.60680    1 45000
## sigma2[8983] 42.38942  8.8424960 28.503292 41.294313  62.75395    1 45996
## sigma2[9021] 28.43124  5.6249579 19.467669 27.698411  41.51889    1 45752
## sigma2[9104] 31.97255  6.3730835 21.822589 31.165836  46.58903    1 46148
## sigma2[9158] 41.60217  8.5227742 28.142369 40.524097  61.46057    1 45000
## sigma2[9198] 33.92014  9.5032647 20.123169 32.280653  56.71477    1 44568
## sigma2[9225] 48.51921 12.2673108 30.178061 46.697395  77.66366    1 45000
## sigma2[9292] 58.08241 21.7665828 29.372380 53.539945 112.85991    1 44575
## sigma2[9340] 51.31195 14.7098521 30.190283 48.792370  86.99259    1 43906
## sigma2[9347] 35.51519  6.9232940 24.439357 34.678935  51.41057    1 45000
## sigma2[9359] 51.88227 10.5635864 35.086644 50.541015  76.39164    1 44555
## sigma2[9397] 44.41285  9.6868769 29.379185 43.050568  66.96302    1 45000
## sigma2[9508] 44.26990 11.4188770 27.390256 42.442255  71.34067    1 45000
## sigma2[9550] 66.29829 18.9323145 38.784889 63.145955 112.07059    1 45717
## sigma2[9586] 42.60858  8.1713380 29.471864 41.628897  61.22517    1 45495
## tau00         9.08284  1.1290112  7.115918  8.999612  11.50486    1 42995
#summary(output)

การเปรียบเทียบโมเดลแบบเบส์สามารถทำได้หลายวิธีการ วิธีการหนึ่งคือการใช้ดัชนี DIC ซึ่ง JAGS คำนวณให้ DIC มีสูตรคำนวณดังนี้

กำหนดให้ deviance มีค่าเท่ากับ\(D(\theta)=-2log(p(\bf{y}|\theta))+C\) โดยที่ \(C\) คือค่าคงที่

การคำนวณค่า DIC มีหลายวิธีการ วิธีการหลักนำเสนอโดย Spiegelhalter และคณะ (2002) และ Gelman และคณะ (2004) ดังนี้

  • Spiegelhalter et al. (2002) กำหนดให้ \(p_D=\overline{D(\theta})-D(\overline{\theta})\)

  • Gelman และคณะ (2004) กำหนดให้ \(p_D=\frac{1}{2}\overline{Var(D(\theta))}\)

\(DIC = p_d+\overline{D(\theta)}\)

หรือ

\(DIC = D(\overline{\theta})+2p_D\)

จากสูตรในข้างต้นจะเห็นว่าโมเดลที่มีค่า DIC ต่ำกว่าจะเป็นโมเดลที่มีความเหมาะสมมากกว่า

dic.samples(fit.null, n.iter=30000, thin = 2)
## Mean deviance:  46741 
## penalty 145.5 
## Penalized deviance: 46887
dic.samples(fit.hetero, n.iter=30000, thin = 2)
## Mean deviance:  46619 
## penalty 320.9 
## Penalized deviance: 46940

package-brms

ผู้เรียนจะเห็นว่าการระบุโมเดลในภาษา JAGS ข้างต้นจำเป็นต้องอาศัยความเชี่ยวชาญพอสมควร นอกจากนี้ยังค่อนข้างใช้เวลาในการเขียนโมเดลถึงแม้จะเป็นโมเดลพื้นฐานอย่าง one-way ANOVA ในข้างต้น ปัจจุบันมี package หลายตัวที่ทำหน้าที่เป็นส่วนต่อประสานระหว่างผู้ใช้กับ engine ต่าง ๆ ซึ่งใช้ภาษาที่เรียบง่ายกว่า ทำให้การ modelling ด้วยวิธีการแบบเบส์ทำได้สะดวกและง่ายขึ้น หัวข้อนี้จะกล่าวถึง package-brms (ย่อมาจาก bayesian multilevel models using stan) ซึ่งเป็น high-level API ตัวหนึ่งของภาษา Stan รายละเอียดมีดังนี้

ขอบเขตของ model ใน package-brms

โมเดลทั่วไปที่สามารถวิเคราะห์ได้ด้วย package-brms มีรูปแบบดังนี้

\(y_i \sim D(f(\eta_i),\theta)\)

เมื่อ \(y_i\) คือค่าสังเกตของตัวแปรตาม ที่มีการแจกแจง \(D\) ในโปรแกรม R จะเรียก \(D\) นี้ว่า family ในขณะที่ \(f(.)\) คือ link function, \(\theta\) คือพารมิเตอร์ของการแจกแจง \(D\) ซึ่งสามารถมีได้หลายตัวขึ้นอยู่กับการแจกแจงที่กำหนด และ \(eta_i\) คือผลรวมเชิงเส้นของตัวแปรอิสระซึ่งสามารถเขียนได้ในรูปทั่วไปดังนี้

\(\bf{\eta}=\bf{X}\beta+\bf{Z}u\)

\(\beta\) คือสัมประสิทธิ์ความถดถอยในระดับ individual ส่วน \(u\) คือสัมประสิทธิ์ความถดถอยในระดับกลุ่ม \(X\) และ \(Z\) คือ design matrix ของตัวแปรอิสระในระดับ individual และ กลุ่ม ตามลำดับ

prior distribution

นอกจากจะมีประสิทธิภาพที่ดีในด้านของอัลกอริทึม MCMC แล้ว Stan ยังสามารถกำหนด prior distribution ได้มีประสิทธิภาพมากกว่าใน JAGS รวมถึง BUGS ด้วย รายละเอียดไปอ่านเอง ไม่อ่านก็ตามใจ

individual-level reg coef —> flat prior

group-level reg coef —> \(\bf{u} \sim N(\bf{0},\Sigma)\)

เมื่อ \(\Sigma\) คือเมทริกซ์ความแปรปรวนร่วม ซึ่งเป็นไปได้ทั้ง diagonal matrix และ symmetry matrix ขึ้นอยู่กับการกำหนด prior distribution ของเมทริกซ์ความแปรปรวนร่วมดังกล่าว นอกจากนี้ยังสามารถ model ให้ group-level reg coef นี้มีการแจกแจงที่เป็นอิสระไปตามกลุ่มได้อีกด้วย —> \(\bf{u}_j \sim N(\bf{0},\Sigma_j)\)

โดยปกติ prior distribution ของ covariance matrix จะกำหนดให้เป็นการแจกแจงแบบ Inverse-Wishart ซึ่งเป็น conjugacy prior กับ exponential family ต่าง ๆ และทำให้อัลกอริทึม Gibb-samplers มีประสิทธิภาพสูง อย่างไรก็ตามใน Stan ไม่จำเป็นต้องกำหนดการแจกแจงในลักษณะดังกล่าว ทำให้การกำหนด prior ของพารามิเตอร์นี้ทำได้ง่ายและมีความหมายที่เข้าใจได้ชัดเจนมากขึ้น ดังนี้

กำหนดให้ \(\Sigma_k=D(\sigma_k)\Omega_kD(\sigma_k)\)

เมื่อ \(D(\sigma_k)\) คือ diagonal matrix ของส่วนเบี่ยงเบนมาตรฐานของ \(u_k\) ส่วน \(\Omega\) คือเมทริกซ์สหสัมพันธ์ของพารามิเตอร์ \(\bf{u}\) ดังกล่าว Lewandowski และคณะ (2009) เสนอให้กำหนด prior ของเมทริกซ์สหสัมพันธ์นี้เป็น LKJ-Correlation prior ที่มีพารามิเตอร์ \(\zeta>0\) กล่าวคือ \(\Omega \sim LKJ(\zeta)\)

  • ถ้า \(\zeta=1\) (ค่าเริ่มต้น) การแจกแจงจะมีลักษณะเป็น uniform บน correlation matrix

  • ถ้า \(\zeta>1\) การแจกแจงจะมีลักษณะลู่เข้าหา identity matrix ขึ้นอยู่กับค่าของ \(zeta\)

  • ถ้า \(0<\zeta<1\) การแจกแจงมีลักษณะเป็น U-shape กล่าวคือให้ค่าความน่าจะเป็นที่จะมีค่า correlation สูง

ส่วน \(sigma_k\) สามารถกำหนด prior distribution แยกรายตัวหรือให้เหมือนกันทั้งหมดก็ได้ ค่าเริ่มต้นของ brms คือการแจกแจงทีแบบครึ่งเดียว (half-Student-t prior) ที่มีองศาความเป็นอิสระเท่ากับ 3

brms ยังสามารถ model ให้พารามิเตอร์ในระดับ individual กับ group มีความสัมพันธ์กันได้ด้วย โดยการแตกเมทริกซ์ \(\Sigma_k\) เป็นดังนี้ \(\Sigma_k=V_k \otimes A_k\) เมื่อ \(V_k\) คือเมทริกซ์ความแปรปรวนร่วมในระดับกลุ่ม (แทน \(\Sigma_k\) ตัวเดิม) และ \(A_k\) คือเมทริกซืความแปรปรวนร่วมของพารามิเตอร์ระหว่างระดับ individual กับ group

การประมาณค่าพารามิเตอร์ในโมเดล

อย่างที่กล่าวไว้ก่อนหน้าแล้วว่า brms ที่ใช้ Stan เป็น engine นั้นมีอัลกอริทึมประมาณค่าพารามิเตอร์ที่มีประสิทธิภาพสูงกว่า JAGS และ BUGS ในแง่ของคุณภาพของตัวอย่างพารามิเตอร์ที่ได้จาก MCMC กล่าวคือตัวอย่างที่ได้จะลู่เข้าหาการแจกแจงความน่าจะเป็นภายหลังได้อย่างรวดเร็ว และมีความค่าอัตสหสัมพันธ์ต่ำ ทำให้ผู้วิเคราะห์ไม่จำเป็นต้องรันลูกโซ่ยาวมากเหมือนกับ Gibb-sampler ใน JAGS และ BUGS อย่างไรก็ตามหากเปรียบเทียบกันต่อรอบ Stan จะใช้เวลามากกว่า

ข้อดีอีกประการหนึ่งของการใช้ Stan ผ่าน brms คือจะให้ค่าสถิติสำหรับเปรียบเทียบโมเดลไว้หลายตัวได้แก่ WAIC และ LOO ซึ่งเป็นดัชนีทีปรับปรุงจาก DIC

ผู้เรียนได้ติดตั้ง package-brms ลงในเครื่องแล้ว ตัวอย่างต่อไปนี้จะแสดงการรัน null model ด้วย brms

library(brms)
fit.brm1<-brm(formula = mathach ~ 1 + (1|schoolid), data=dat1,
              family=gaussian(),
              warmup=1000,
              iter=3000,
              chains=3)
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## 
## SAMPLING FOR MODEL '0ead22ac663f58d23a13a947854bc20f' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000481 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 4.81 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 3000 [  0%]  (Warmup)
## Chain 1: Iteration:  300 / 3000 [ 10%]  (Warmup)
## Chain 1: Iteration:  600 / 3000 [ 20%]  (Warmup)
## Chain 1: Iteration:  900 / 3000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1001 / 3000 [ 33%]  (Sampling)
## Chain 1: Iteration: 1300 / 3000 [ 43%]  (Sampling)
## Chain 1: Iteration: 1600 / 3000 [ 53%]  (Sampling)
## Chain 1: Iteration: 1900 / 3000 [ 63%]  (Sampling)
## Chain 1: Iteration: 2200 / 3000 [ 73%]  (Sampling)
## Chain 1: Iteration: 2500 / 3000 [ 83%]  (Sampling)
## Chain 1: Iteration: 2800 / 3000 [ 93%]  (Sampling)
## Chain 1: Iteration: 3000 / 3000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 9.68877 seconds (Warm-up)
## Chain 1:                10.1136 seconds (Sampling)
## Chain 1:                19.8024 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '0ead22ac663f58d23a13a947854bc20f' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000322 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 3.22 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 3000 [  0%]  (Warmup)
## Chain 2: Iteration:  300 / 3000 [ 10%]  (Warmup)
## Chain 2: Iteration:  600 / 3000 [ 20%]  (Warmup)
## Chain 2: Iteration:  900 / 3000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1001 / 3000 [ 33%]  (Sampling)
## Chain 2: Iteration: 1300 / 3000 [ 43%]  (Sampling)
## Chain 2: Iteration: 1600 / 3000 [ 53%]  (Sampling)
## Chain 2: Iteration: 1900 / 3000 [ 63%]  (Sampling)
## Chain 2: Iteration: 2200 / 3000 [ 73%]  (Sampling)
## Chain 2: Iteration: 2500 / 3000 [ 83%]  (Sampling)
## Chain 2: Iteration: 2800 / 3000 [ 93%]  (Sampling)
## Chain 2: Iteration: 3000 / 3000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 9.70008 seconds (Warm-up)
## Chain 2:                10.0076 seconds (Sampling)
## Chain 2:                19.7077 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '0ead22ac663f58d23a13a947854bc20f' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000372 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 3.72 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 3000 [  0%]  (Warmup)
## Chain 3: Iteration:  300 / 3000 [ 10%]  (Warmup)
## Chain 3: Iteration:  600 / 3000 [ 20%]  (Warmup)
## Chain 3: Iteration:  900 / 3000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1001 / 3000 [ 33%]  (Sampling)
## Chain 3: Iteration: 1300 / 3000 [ 43%]  (Sampling)
## Chain 3: Iteration: 1600 / 3000 [ 53%]  (Sampling)
## Chain 3: Iteration: 1900 / 3000 [ 63%]  (Sampling)
## Chain 3: Iteration: 2200 / 3000 [ 73%]  (Sampling)
## Chain 3: Iteration: 2500 / 3000 [ 83%]  (Sampling)
## Chain 3: Iteration: 2800 / 3000 [ 93%]  (Sampling)
## Chain 3: Iteration: 3000 / 3000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 9.15719 seconds (Warm-up)
## Chain 3:                9.2879 seconds (Sampling)
## Chain 3:                18.4451 seconds (Total)
## Chain 3:
plot(fit.brm1)

summary(fit.brm1)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: mathach ~ 1 + (1 | schoolid) 
##    Data: dat1 (Number of observations: 7185) 
##   Draws: 3 chains, each with iter = 3000; warmup = 1000; thin = 1;
##          total post-warmup draws = 6000
## 
## Group-Level Effects: 
## ~schoolid (Number of levels: 160) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     2.95      0.19     2.60     3.33 1.00     1294     2477
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    12.64      0.25    12.17    13.15 1.00      705     1210
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     6.26      0.05     6.15     6.37 1.00    12329     4130
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

posterior predictive check

pp_check(fit.brm1, resp = "mathach")
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

bayes_R2(fit.brm1)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.1725417 0.007560215 0.1576866 0.1873903
loo(fit.brm1)
## 
## Computed from 6000 by 7185 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo -23445.7 49.4
## p_loo       146.3  2.1
## looic     46891.4 98.9
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

เราสามารถปรับเปลี่ยน prior ได้หลายวิธีการ วิธีการหนึ่งคือการเรียกค่าเริ่มต้นของ prior มาก่อน แล้วปรับเปลี่ยนทีละตัวตามต้องการดังนี้

priors<-get_prior(mathach ~ 1 + (1|schoolid), data=dat1,
              family=gaussian())
priors
##                    prior     class      coef    group resp dpar nlpar bound
##  student_t(3, 13.1, 8.1) Intercept                                         
##     student_t(3, 0, 8.1)        sd                                         
##     student_t(3, 0, 8.1)        sd           schoolid                      
##     student_t(3, 0, 8.1)        sd Intercept schoolid                      
##     student_t(3, 0, 8.1)     sigma                                         
##        source
##       default
##       default
##  (vectorized)
##  (vectorized)
##       default
priors$prior[1]<-"normal(0,20)" #mu and sd

จากนั้นก็ run ใหม่โดยกำหนดอาร์กิวเมนท์ prior ดังนี้

fit.brm2<-brm(formula = mathach ~ 1 + (1|schoolid), data=dat1,
              family=gaussian(),
              prior=priors,
              warmup=1000,
              iter=3000,
              chains=3)
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## 
## SAMPLING FOR MODEL '2b3230ad39a0e3f7602f99f6d472f68b' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000464 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 4.64 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 3000 [  0%]  (Warmup)
## Chain 1: Iteration:  300 / 3000 [ 10%]  (Warmup)
## Chain 1: Iteration:  600 / 3000 [ 20%]  (Warmup)
## Chain 1: Iteration:  900 / 3000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1001 / 3000 [ 33%]  (Sampling)
## Chain 1: Iteration: 1300 / 3000 [ 43%]  (Sampling)
## Chain 1: Iteration: 1600 / 3000 [ 53%]  (Sampling)
## Chain 1: Iteration: 1900 / 3000 [ 63%]  (Sampling)
## Chain 1: Iteration: 2200 / 3000 [ 73%]  (Sampling)
## Chain 1: Iteration: 2500 / 3000 [ 83%]  (Sampling)
## Chain 1: Iteration: 2800 / 3000 [ 93%]  (Sampling)
## Chain 1: Iteration: 3000 / 3000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 8.98691 seconds (Warm-up)
## Chain 1:                9.02464 seconds (Sampling)
## Chain 1:                18.0116 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '2b3230ad39a0e3f7602f99f6d472f68b' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000285 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 2.85 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 3000 [  0%]  (Warmup)
## Chain 2: Iteration:  300 / 3000 [ 10%]  (Warmup)
## Chain 2: Iteration:  600 / 3000 [ 20%]  (Warmup)
## Chain 2: Iteration:  900 / 3000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1001 / 3000 [ 33%]  (Sampling)
## Chain 2: Iteration: 1300 / 3000 [ 43%]  (Sampling)
## Chain 2: Iteration: 1600 / 3000 [ 53%]  (Sampling)
## Chain 2: Iteration: 1900 / 3000 [ 63%]  (Sampling)
## Chain 2: Iteration: 2200 / 3000 [ 73%]  (Sampling)
## Chain 2: Iteration: 2500 / 3000 [ 83%]  (Sampling)
## Chain 2: Iteration: 2800 / 3000 [ 93%]  (Sampling)
## Chain 2: Iteration: 3000 / 3000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 9.12027 seconds (Warm-up)
## Chain 2:                8.97629 seconds (Sampling)
## Chain 2:                18.0966 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '2b3230ad39a0e3f7602f99f6d472f68b' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000284 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 2.84 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 3000 [  0%]  (Warmup)
## Chain 3: Iteration:  300 / 3000 [ 10%]  (Warmup)
## Chain 3: Iteration:  600 / 3000 [ 20%]  (Warmup)
## Chain 3: Iteration:  900 / 3000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1001 / 3000 [ 33%]  (Sampling)
## Chain 3: Iteration: 1300 / 3000 [ 43%]  (Sampling)
## Chain 3: Iteration: 1600 / 3000 [ 53%]  (Sampling)
## Chain 3: Iteration: 1900 / 3000 [ 63%]  (Sampling)
## Chain 3: Iteration: 2200 / 3000 [ 73%]  (Sampling)
## Chain 3: Iteration: 2500 / 3000 [ 83%]  (Sampling)
## Chain 3: Iteration: 2800 / 3000 [ 93%]  (Sampling)
## Chain 3: Iteration: 3000 / 3000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 9.80854 seconds (Warm-up)
## Chain 3:                9.52466 seconds (Sampling)
## Chain 3:                19.3332 seconds (Total)
## Chain 3:

ผู้วิเคราะห์สามารถทดสอบสมมุติฐานของพารามิเตอร์ที่ต้องการในโมเดลได้ โดยใช้วิธีการเปรียบเทียบ test value กับ HDI โดยใช้ฟังก์ชัน hypothesis() ดังตัวอย่างต่อไปนี้

hypothesis(fit.brm1, "Intercept>0", class="sd", group="schoolid", alpha=0.05)
## Hypothesis Tests for class sd_schoolid:
##        Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (Intercept) > 0     2.95      0.19     2.65     3.27        Inf         1
##   Star
## 1    *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Extract output from brm

ผู้วิเคราะห์สามารถแปลงผลการวิเคราะห์ที่ได้ใช้เป็น mcmc objects แล้วไปวิเคราะห์แบบเดิมก็ได้ดังนี้

gamma.mcmc<-as.mcmc(as.data.frame(fit.brm1, pars="b_Intercept"))
## Warning: Argument 'pars' is deprecated. Please use 'variable' instead.
tau.mcmc<-as.mcmc(as.data.frame(fit.brm1, pars="sd")^2)
## Warning: Argument 'pars' is deprecated. Please use 'variable' instead.
sigma2.mcmc<-as.mcmc(as.data.frame(fit.brm1, variable="sigma"))^2
icc<-tau.mcmc/(tau.mcmc+sigma2.mcmc)
mean(icc)
## [1] 0.1820599
sd(icc)
## [1] 0.01892475
plotPost(icc, compVal=0.05, cenTend="median", ROPE=c(0.04,0.06), credMass=0.95)

##                  ESS      mean    median     mode hdiMass    hdiLow   hdiHigh
## Param. Val. 1323.766 0.1820599 0.1814084 0.181766    0.95 0.1445473 0.2178896
##             compVal pGtCompVal ROPElow ROPEhigh pLtROPE pInROPE pGtROPE
## Param. Val.    0.05          1    0.04     0.06       0       0       1

อีกลักษณะหนึ่งสามารถทำได้โดยใช้ package-tidybayes ช่วย สามารถอ่านเพิ่มเติมได้จาก ไม่อ่านก็ตามใจ

library(tidybayes)
## 
## Attaching package: 'tidybayes'
## The following objects are masked from 'package:brms':
## 
##     dstudent_t, pstudent_t, qstudent_t, rstudent_t
get_variables(fit.brm1)%>%head(15)
##  [1] "b_Intercept"                "sd_schoolid__Intercept"    
##  [3] "sigma"                      "r_schoolid[1224,Intercept]"
##  [5] "r_schoolid[1288,Intercept]" "r_schoolid[1296,Intercept]"
##  [7] "r_schoolid[1308,Intercept]" "r_schoolid[1317,Intercept]"
##  [9] "r_schoolid[1358,Intercept]" "r_schoolid[1374,Intercept]"
## [11] "r_schoolid[1433,Intercept]" "r_schoolid[1436,Intercept]"
## [13] "r_schoolid[1461,Intercept]" "r_schoolid[1462,Intercept]"
## [15] "r_schoolid[1477,Intercept]"
fit.brm1%>%spread_draws(r_schoolid[schoolid, Intercept]) #group-mean
## # A tibble: 960,000 × 6
## # Groups:   schoolid, Intercept [160]
##    schoolid Intercept r_schoolid .chain .iteration .draw
##       <int> <chr>          <dbl>  <int>      <int> <int>
##  1     1224 Intercept      -2.09      1          1     1
##  2     1224 Intercept      -2.74      1          2     2
##  3     1224 Intercept      -3.97      1          3     3
##  4     1224 Intercept      -1.75      1          4     4
##  5     1224 Intercept      -3.37      1          5     5
##  6     1224 Intercept      -1.55      1          6     6
##  7     1224 Intercept      -2.91      1          7     7
##  8     1224 Intercept      -3.87      1          8     8
##  9     1224 Intercept      -1.81      1          9     9
## 10     1224 Intercept      -2.17      1         10    10
## # … with 959,990 more rows
fit.brm1%>%spread_draws(b_Intercept,sd_schoolid__Intercept) #grand-mean and sqrt(tau00)
## # A tibble: 6,000 × 5
##    .chain .iteration .draw b_Intercept sd_schoolid__Intercept
##     <int>      <int> <int>       <dbl>                  <dbl>
##  1      1          1     1        12.5                   2.98
##  2      1          2     2        12.3                   3.01
##  3      1          3     3        12.5                   3.14
##  4      1          4     4        12.5                   3.00
##  5      1          5     5        12.6                   2.92
##  6      1          6     6        12.8                   2.72
##  7      1          7     7        12.5                   2.94
##  8      1          8     8        12.5                   3.15
##  9      1          9     9        12.1                   2.90
## 10      1         10    10        12.3                   2.76
## # … with 5,990 more rows
# mean and median with quantile interval
fit.brm1%>%spread_draws(b_Intercept,sd_schoolid__Intercept)%>%
              mean_qi()
## # A tibble: 1 × 9
##   b_Intercept b_Intercept.lower b_Intercept.upper sd_schoolid__Intercept
##         <dbl>             <dbl>             <dbl>                  <dbl>
## 1        12.6              12.2              13.1                   2.95
## # … with 5 more variables: sd_schoolid__Intercept.lower <dbl>,
## #   sd_schoolid__Intercept.upper <dbl>, .width <dbl>, .point <chr>,
## #   .interval <chr>
fit.brm1%>%spread_draws(b_Intercept,sd_schoolid__Intercept)%>%
              median_qi()
## # A tibble: 1 × 9
##   b_Intercept b_Intercept.lower b_Intercept.upper sd_schoolid__Intercept
##         <dbl>             <dbl>             <dbl>                  <dbl>
## 1        12.6              12.2              13.1                   2.94
## # … with 5 more variables: sd_schoolid__Intercept.lower <dbl>,
## #   sd_schoolid__Intercept.upper <dbl>, .width <dbl>, .point <chr>,
## #   .interval <chr>

random-coefficient model

priors<-get_prior(mathach~1+ses+(1+ses|schoolid), 
            data=dat, family=gaussian())
priors
##                    prior     class      coef    group resp dpar nlpar bound
##                   (flat)         b                                         
##                   (flat)         b       ses                               
##                   lkj(1)       cor                                         
##                   lkj(1)       cor           schoolid                      
##  student_t(3, 13.1, 8.1) Intercept                                         
##     student_t(3, 0, 8.1)        sd                                         
##     student_t(3, 0, 8.1)        sd           schoolid                      
##     student_t(3, 0, 8.1)        sd Intercept schoolid                      
##     student_t(3, 0, 8.1)        sd       ses schoolid                      
##     student_t(3, 0, 8.1)     sigma                                         
##        source
##       default
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
priors$prior[1]<-"normal(0,20)" #mu and sd

fit.brm3<-brm(mathach~1+ses+(1+ses|schoolid), 
            data=dat, family="gaussian",
            prior=c(set_prior("horseshoe(1)", class="b"),
                    set_prior("student_t(1,0,5)", class="sd"),
                    set_prior("lkj(0.8)", class="cor")),
            warmup=1000, iter=6000, chains=3, thin=5,
            cores=3)
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Warning: There were 323 divergent transitions after warmup. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: The largest R-hat is 1.06, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#r-hat
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#tail-ess
plot(fit.brm3)

summary(fit.brm3)
## Warning: There were 323 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See http://mc-stan.org/misc/
## warnings.html#divergent-transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: mathach ~ 1 + ses + (1 + ses | schoolid) 
##    Data: dat (Number of observations: 7185) 
##   Draws: 3 chains, each with iter = 1200; warmup = 200; thin = 5;
##          total post-warmup draws = 600
## 
## Group-Level Effects: 
## ~schoolid (Number of levels: 160) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)          2.21      0.15     1.94     2.53 1.00      974     2758
## sd(ses)                0.53      0.23     0.04     0.91 1.00      612     1852
## cor(Intercept,ses)    -0.16      0.29    -0.83     0.41 1.00     1929     2009
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    12.63      0.21    12.25    13.02 1.04       49       79
## ses           2.38      0.12     2.15     2.62 1.00     1471     2201
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     6.07      0.05     5.97     6.17 1.01      238      190
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
pp_check(fit.brm3, resp = "mathach")

bayes_R2(fit.brm3)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2206894 0.008245877 0.2042388 0.2371053
loo(fit.brm3)
## 
## Computed from 3000 by 7185 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo -23238.5  51.2
## p_loo       160.3   2.5
## looic     46476.9 102.4
## ------
## Monte Carlo SE of elpd_loo is 1.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.